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I.  INTRODUCTION 


Vortices  and  vortex  wakes  have  become  a  major  theme  of  aerodynamics  research 
since  the  advent  of  the  large  aircraft  and  the  understanding  of  their  evolution  required 
an  examination  of  many  fundamental  problems  in  .fluid  mechanics.  Much  of  the 
progress  made  during  the  past  two  decades  was  discussed  at  the  Symposium  on 
Aircraft  Wake  Turbulence  and  Its  Detection  [Ref  1]  and  at  the  Aircraft  Wake  Vortices 
Conference  [Ref  2].  Comprehensive  reviews  of  the  entire  subject  have  been  given  by 
Donaldson  and  Bilanin  [Ref  3],  Widnall  [Ref  4]  and  Hallock  and  Eberle  [Ref  5]. 

These  studies,  as  well  as  numerous  others  carried  out  since  1977,  have  uncovered 
a  number  of  complex  problems  which  must  be  resolved  in  order  to  achieve  a  better 
understanding  of  the  important  features  of -trailing  vortices  in  homogeneous  and 
stratified  media.  The  principal  ones  are  as  follows  [Ref  6]: 

(a)  Roll-up  process;  The  velocity  and  turbulence  distribution  at  anv  station 
behind  the  wing  depend  on  the  wing  section,  wins-tip  shape^.  Revnolds 
number,  wins  incidence,  and  the  distance  of  the  stationTrom  the  wing  [Ref  7]. 
The  distribufions  of  the  initial  velocitv  and  turbulence,  which  influence  the 
roll-up  and  the  decav  process,  cannot  be  changed  independently.  For 
example,  a  chanse  in  tip  shape  changes  the  core  size,  as  well  .as  the  velocity 
and  turbulence  distributions.  Hish  levels  of  turbulence  result  m  an  increased 
diffusion  of  vorticity,  which  in  turn  increases  the  core  size. 

(b)  Probe  sensitivity  of  the  vortices:  Flow  visualization  studies  suggest  that 
trailing  vortices 'are  extremelv  sensitive  to  disturbances  created  .by  "even  veiw 
small  probes  or  bubbles.  This  forces  one  to  use  non-intrusive  means .  of 
measurements  such  as  an  LDV.  Even  then,  'vortex  wandering'  [Ref  8],  which 
makes  the  vortices  appear  larger  than  normal  in  time-averaged  velocity 
measurements  (for  vortices  senerated  by  a  wing  in  a  wind  tunnel),  or  the 
unsteadv  nature  of  the  flow  ffor  vortices' senerated  bv  a  wing  in  a  tow  basin) 
makes  the  mean  velocity  prohles  in  the  vortices  difficult  to  determine. 

(c)  Large-scale  instabilities:  The  vortices  are  seldom  observed  to  decay  aw-ay 
•  owihs  to  viscous  and  turbulent  dissipation,  but  are  almost  always. 'destroyed 

bv  either  mutual  induction  instability  (Crow  instability  [Ref  9])  and/or  vortex 
breakdown.  The  Crow  instability  grows  exponentiallv,  and  results  either  in  a 
linking  of  the  vortex  pair  into  a  series  of  crude  vortex  rings  or  in  a  hishly 
disorsanized  interminslins  of  the  vortices.  Vortex  breakdown,  wlios'e 
mathematical  details  have"  not  yet  been  adeouatelv  treated,  rearranges  the 
vortex  structure  and  increases  the'  core  size,  turbulence,  and  energy  dissipation. 
Thus,  it  is  verv  difficult  to  rneasure  accurately  th.e  trajector.ies"bf  the  three- 
dimensional  vo'rtices  from  their  creation  to  their  ultimate  demise. 

(d)  Reynolds  number;  Even  the  highest  Revnolds  numbers,  based  on  wimi  chord, 
rea'ched  in  wind  tunnels  or  tovinsz  basins,  are  in  order  of  maanatude  lower 
than  what  is  possible  for  an  aircraft.  Thus,  the  scale  effects  are  not  easy  to 
assess. 

(e)  Ambient  conditions  such  as  turbulence  and  stratification  play  major  roles  in 
the  evolution  of  vortices.  The  quantification  of  these  'effects  requires 
numerical  analysis  and  extremely  careful  experiments. 
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(f)  Ground  or  free-surface  efTects:  The  vortex,  pair  mav  move  towards  a  risid 
boundary'  at  which,  the  no-slip  condition  must  be  satisfied  or  towards  a  free 
surface  at  which  the  zero-shear  condition  must  be  satisfied.  In  either  case,  the 
vortices  come  under  the  influence  of  their  images  and  move  accordingly. 

The  phenomenon  is  further  complicated  by  several  additional  facts.  When  the 
vortices  are  propelled  towards  a  rigid  surface,  vorticity  of  opposite  sign  is  generated  on 
the  no-slip  boundary*  and  swept  towards  the  vortex  pair.  The  total  vorticity  diminishes 
very  quickly  as  vorticity  from  the  two  regions  diffuses,  the  wall  region  serving  as  a 
strong  sink  for  the  vorticity  associated  with  the  -original  vortex  [Ref  10].  The 
developement  of  a  boundary  layer  along  the  rigid  wall  may  give  rise  to  flow  separation 
for  sufficiently  high  Reynolds  numbers.  With  or  without  such  a  separation,  however, 
the  center  of  the  vortex  pair  eventually  moves  away,  or' 'rebounds',  from  the  wall 
[Refs.  10-12]. 

For  the  case  of  a  zero-stress  boundary,  the  free  surface  still  acts  as  a  vorticity 
sink,  but  this  is  relatively  weak  due  to  the  absence  of  intense  oppositely-signed 
vorticity.  Thus,  in  the  absence  of  other  impeding  phenomena,  one  expects  a  mild 
interaction  betw'een  the  vortices  and  the  free  surface  and  a  small  rebound  of  the  vortex 
center  from  the  free  surface.  However,  the  ability  of  the  free  surface  to  deform  under' 
the  influence  of  strain  fields  leads  to  a  strong  interaction  between  the  vortices  and  the 
free  surface. 

It  is  evident  from  the  foregoing  that  the  motion  and  the  life-span  of  trailing, 
vortices  are  goverened  by  a  number  of  nonlinearly-dependent  complex  phenomena.'  A 
number  of  experimental  and  analytical  studies  have  been  carried  out  at  the  Naval 
Postgraduate  School  by  Sarpkaya  and  his  students  [Refs.  6,13-17]  in  order  to 
investigate  the  efTects  of  these  parameters  on  the  rise  and  demise  of  the  trailing  vortices 
in  homogeneous  and  density  stratified  media.  These  studies  have  clearly  identified  the 
various  demise  mechanisms  in  both  media  and  established  basic  relationships  between 
.the  rise  of  vortices  and  the  governing  parameters  in  a  finite  as  well  as  effectively  infinite 
medium  [Ref  IS],  free  from  ambient  turbulence.  The  present  investigation  is  a 
continuation  of  the  foregoing  studies  and  is  confined  to  the  computer  modeling,  via  a 
finite  differencing  technique,  of  the  rise  and  demise  of  trailing  vortices  in  stratified 
media. 
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11.  EXPERIMENTAL  EQUIPMENT  AND  PROCEDURES 


A.  EQUIPMENT 

The  equipment  used  to  generate  the  trailing  vortices  has  been  extensively  used  at 
this  facility  over  the  past  three  years  [Refs.  6,14-19].  Only  the  salient  features,  most 
recent  modifications  and  the  adaptation  for  this  work  are  briefly  described  in  the 
following. 

The  system  consists  essentially  of  a  towing  basin.  The  auxiliary  components 
consist  of  the  brine  fill  system,  turbulence  management  system,  top  and  bottom 
carriages,  velocity  measuring  system,  hghting  system,  and  the  model  [Refs.  14-19].  The 
brine  fill  system  is  controlled  by  a  computer  to  provide  a  consistently  repeatable 
density  stratified  media  [Ref.  19]. 

Drains  are  provided  at  the  bottom  of  each  end  of  the  basin.  Two  parallel  rails 
are  mounted  along  the  bottom  of  the  tank.  A  carriage  rides  smoothly  on  these  rails 
and  provides  the  test  body  with  a  constant  velocity  through  the  use  of  an  endless  cable 
and  a  DC  motor.  The  velocity  of  the  model  is  measured  and  continuously  monitored 
through'  the  use  of  a  magnetic  linear  displacement  transducer.  It  yields  a  signal 
proportional  to  the  velocity  of  the  model  with  an  error  of  less  than  one  percent. 

The  two  rails,  the  carriage  and  the  filling  pipes  are  located  on  or  near  the  bottom 
of  the  basin  and  under  a  turbulence  management  system.  This  system  consists  of  one 
inch  thick  polyurethane  foam  sandwiched  between  two  perforated  aluminum  plates. 

B.  MODEL 

.A  rectangular  foil  (NACA  0012)  was  used  in  the  present  study.  The  fojl  had  the 
following  dimensions:  B  =  4.50  in.,  b^  =  3.53  in.,  and  c  =  2.32  in.  The  interior  of 
this  model  vvas  hollowed  and  used  as  dye  reservoir  to  seed  the  vortex  cores.  The  model 
was  mounted  on  its  base  by  means  of  a  thin  streamlined  aluminum  bar  with  a  cross 
section  of  a  NACA  0006  foil  and  set  at  the  desired  angle  of  attack.  As  noted  earlier, 
the  model  was  pulled  by  means  of  a  DC  motor,  pulley,  and  cable  system  at  the  desired 
speed.  Additional  details  of  the  model  construction  and  mounting  are  discussed  by 
Johnson  [Ref.  15]. 


C  TEST  PROCEDURES 


A  model  is  placed  in  the  basin  and  the  basin  is  filled  gradually  with  a  brine 
solution  [Ref.  19]  to  the  desired  le\-el.  In  the  cases  where  dye  was  used  to  seed  the 
vortex  pair,  the  model  was  filled  with  dye  prior  to  filling  the  basin.  .After  removal  of 
trapped  air  and  after  sufficient  time  to  allow  for  the  escape  of  dissolved  air  in  the  water 
and  the  elimination  of  any  internal  currents  in  the  basin,  the  model  was  set  in  motion. 
Most  of  the  experiments  were  repeated  at  least  three  times. 

The  motion  of  the  trailing  vortices  were  recorded  on  highspeed  Polaroid  film  at 
the  test  section  (one  of  the  plexiglass  panels  near  the  middle  of  the  basin).  Each 
picture  included  tw'o  clocks  accurate  to  0.1  seconds,  the  vertical  and  horizontal  scales 
on  the  plexiglass  window  and,  of  course,  the  side  view  of  the  trailing  vortices  as  they 
rose  from  the  model  after  formation.  The  time  interval  between  successive  pictures  is 
determined  from  the  two  clocks.  The  first  picture  always  captured  the  model  the 
instant  it  passed  through  the  test  section.  The  subsequent  pictures  (taken  at  about 
0.75  second  intervals)  captured  the  rise  and  demise  of  the  vortices.  The  vertical  rise  of 
the  vortices  is  determined  from  the  verticle  scale.  Attention  has  been  paid  to  the  fact 
that  the  vortices  are  farther  away  from  the  scale  on  the  window  and  that  the  scale 
placed  vertically  in  the  middle  of  the  test  section  does  not  exactly  correspond  to  the 
scale  marked  on  the  window  due  to  refraction  and  parallax.  The  necessary  correction 
was  made  by  photographing  a  scale  placed  in  "water  in  the  middle  of  the  test  section 
together  with  a  scale  marked  on  the  window*.  This  resulted  in  a  simple  conversion 
table  which  enabled  the  determination  of  the  actual  position  of  the  vortex  from  the 
scale  reading  on  the  photograph. 

The  results  are  normalized  and  plotted  in  various  forms  and  compared  with  those 
obtained  in  the  previous  runs.  Each  experiment  w'as  repeated  three  times  for  most  of 
the  data  presented  herein.  All  trailing  vortices  w'ere  recorded  on  film  until  the  time 
they  have  completely  dissipated  either  due  to  aging  (diffusion  of  vorticity  due  to 
viscousity,  turbulence,  entrainment,  and  detrainment),  or  due  to  Crow  instability 
(sinusoidal  instability  and  linking  of  vortices),  and,  or  due  to  vortex  breakdown  (core 
bursting).  It  was  thus"  possible  to  determine  the  life  span  of  vortex  cores  from  their 
formation  to  their  final  demise. 
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D.  COMPUTER  MODEL 

The  computer  model  (Appendix  A)  is  run  on  the  IBM  3300  main-frame  computer 
at  the  Naval  Postgraduate  School.  The  computer  code  was  written  to  run  on  either  the 
MVS  Batch  Processing  system  or  the  VM,CMS  system.  The  stratification  parameter 
specified  for  the  experimental  runs  is  used  to  determine  other  parameters  [Ref  1']  used 
by  the  program.  The  code  is  compiled  under  VS  Fortran  and  is  loaded  and  executed 
using  DISSPLA  version  10.0.  The  code  was  \^'ritten  such  that  the  user  has  the  capacity 
to  stipulate  the  number  of  time  step  iterations  that  are  executed.  .Also,  the  time  steps 
at  which  a  constant  density  contour,  complex  velocity  profile,  density  perturbation 
contour  and  the  vorticity  contour  are  plotted,-  are  designated  by  the  user  prior  to 
compilation.  These  plots  are  analized  to  determine  the  normalized  rise  height  of  the 
initial  vortex  pair  until  its  demise.  These  are  then  compared  with  those  obtained 
experimentally. 

The  computer  model  is  an  intense  computational  algorithm  that  requires 
approxiamately  12.5  minutes  of  Central  Processing  Unit  (CPU)  time  per  iteration.  Due 
to  this  high  use  of  CPU  time,  the  code  has  been  designed  to  he  executed  until  stoped  at 
a  user  defined  step  iteration  and  then  save  all  necessary  data  required  to  restart  the 
model.  The  code  can  then  be  restarted  at  that  step  iteration  at  a  later  time.  The 
program  VORRS  (VORtex  ReStart)  has  been  specifically  w'ritten  to  accomplish  this 
task. 
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III.  ANALYSIS 


A.  NUMERICAL  ANALYSIS  OF  THE  VORTEX  .MOTIONS 

The  seneration  of  the  internal  waves  and  the  rise  and  demise  of  the  vortices  in  a 
stratified  medium  may  well  be  analyzed  through  the  use  of  the  equations  of  motion  for 
an  incompressible  fluid  for  both  laminar  and  turbulent  motions  provided  that  a  suitable 
turbulence  closure  model  is  adopted  and  the  usual  Boussinesq  approximation 
(gravitational  acceleration  is  much  larger  than  the  fluid  accelerations)  is  made.  For  the 
type  of  motions  considered  herein  the  Boussinesq  approximation  is  quite  valid  and  has 
been  used  in  the  investigation  of  all  types  of  internal  waves  in  stratified  fluids. 

For  a  two-dimensional  flow,  with  y  vertical  and  x  horizontal,  the  equations  of 
motion  are 


5u  du  5u  1  5p  •  T 

—  +  u  —  +  V  —  =  - - +  vV^u 

5t  dx  dv  p  dx 


and 


dv  d\  dv  1  5p  , 

- h  u  —  +  V  —  =  -  g  - - +  vV'^v 

5t  dx  dy  p  dy 


(3.1) 


(3.2) 


and  the  equation  of  continuity 


5u  dv 

—  +  -  =  0. 

dx  dy 


(3.3) 


Defining  vorticity  as  usual  by 


3u  dv 
dy  dx 


(3.4) 


and  eliminating  pressure  between  Eq.  3.1  and  Eq.  3.2  ,  one  has 
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(3.5) 


_L  H - Z  H - Z 

di  dx  dy 


g£p' 

Po 


in  which  p'  is  the  fluctuating  part  of  the  density  p  given  by 


P  =  Po  P(y)  ^  P'(x,y,t) 


(3.6) 


where  p^  is  the  reference  density,  and  p(y)  the  initial  density  at  elevation  y  at  t  =  0. 
The  diffusion  of  density  is  given  by 


d  d  d 

( —  +  u  —  +  V  — 
5t  dx  dy 


-  vV^)p  =  0 


(3.7) 


Combining  Eqs.  3.3.  3.6  and  3.7,  and  simplifying,  one  has 


dp'  dup'  dyp' 
—  ^ - -  4-  - - 

di  dx  dy 


V 


—  +  v7-p'  +  V' 
dy 


d^p 



t7y“ 


(3.8) 


It  is  convenient  to  cast  the  foregoing  equations  into  nondimensional  forms, 
scaling  each  variable  by  a  quantity  characteristic  of  its  expected  magnitude.  Two 
possible  time  scales  exist.  There  is  the  dynamic  time  scale,  which  is  the  time  a 
characteristic  length  would  be  traversed  by  a  fluid  particle  traveling  at  the  characteristic 
velocity,  and  there  is  the  buoyant  time  scale  based  on  the  natural  buoyancy  frequency 
of  the  stratified  flow,  i.e.,  the  Brunt-Vaisala  frequency  N,  as  defined  by  Daly  [Ref  17]. 
Each  of  the  scales  gives  a  slightly  different  form  of  the  normalized  qoverning  equations. 

Dynamic  Scale: 

Introducing  and  as  the  characteristic  velocity  and  length,  one  has 


^L,  U  t  U  L  U 

Cm  =  - tm  =  — ,  Re  =  F  =  — 

L  ’  V  ’  (eL)' 


Normalizing  Eqs.  3.5  and  3.8 


m 


5Cm  ^^mCm  5V  C  1  .  1  C*  P 

-  -m  ^  (  m  -m  ^  m  -m  \  _  _  ^2^  j _ _ 

5t  ^Xm  Re  F  -  dx 

ni  m  ^  rn  v  m 


(3.9) 


(3.10) 
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and 


dp  '  p_'  dvo'  ^  1  ,  d^p^ 

Jim  ^  ^ - m£m  j - ^ - ^  -—S'-) 


m 


dx 


m 


0\' 


m 


Re 


m 


OX 


m 


(3.11) 


in  which 


^  U  ’ 

F. .-  =  — ^.  Re  =  -^.  n  =  ^ 


sL. 


L>. 

V 


N 

_ ( 

N 


a 

o 


N/  =  — ,  N  2  = 

C  T  ’  O 


g  p(y) 
p.  ■ 


g  ^Pm  _ 


L  dy 


m 


Equations  3.10  ans  3.11  are  valid  when  >  >  1,  i.e.,  when  the  bouyancy  has 
ver>'  little  influence  on  the  nonlinear  dynamics  of  the  motion.  In  this  case,  the  density 
perturbation  acts  as  a  passive  scalar  advected  by  the  velocity  field. 

Bouyant  scale: 

Introducing  the  following  dimensionless  parameters; 


V  = 
'  m 


Pm  -  ’  TT  *  TT  ' 


m 


= 


u 


m 


V 

u" 


'^'c.  ‘Sm 


C  ■■ 


m 

T  2 
c 


X 

r 


L. 


(3.12) 


and  normalizing  Eqs.  3.5  and  3.8,  one  has 


dC 


du  C  dy 

—  —  m 


~m  p  /  ni  •‘m  m  \ 

dx_  'dy  ’ 


=  h  +  fPni 


in 


m 


m 


Re 


<2X 


m 


(3.13) 


and 


dp^'  ^u_p_  dy  p 


<2t 


+  Fv(- 


m 


dx 


ni 


+  — mi-iii)  =  F  n-v 

dy 


m 


m 


(3.14) 


in  which  F^,,  Re  and  n"  are  as  defined  in  Eq.  3. 1 1. 
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Equations  3.13  and  3.14  are  valid  when  <  <  1  and  buoyancy  dominates  the 

flow.  When  approaches  zero,,  the  equations  that  result  from  Eqs-3.13  and  3.14 

describe  the  propogation  of  linear  internal  waves.  The  F^.  <  <  1  regime  is  of  interest 
in  the  present  investigation  because  for  submerged  bodies  of  naval  interest  F^.  (the 
Froude  number)  is  about  0.01.  The  analysis  will  consider  the  full  nonlinear  equations 
given  by  Eqs.  3.13  and  3.14  rather  than  their  linearized  form  (i.e..  the  case  of  F^.  =  0). 
However,  the  flow  will  be  assumed  to  be  inviscid,  i.e.,  the  viscous  difiusion  will  be 
ignored  to  a  first  order  approximation.  Subsequent  analysis  will  incorporate  the  eflects 
of  viscous  and  turbulent  diffusion  into  the  numerical  calculations.  The  velocities  u  and 
V  are  given  by  the  Biot-Savart  law  as 


u(x,y) 


J 


(y^  -  yK(x\y')dx'dy^ 
27tr^ 


(3.15) 


and 


u(x,y)  =  J 


.(x  -  x')C('x',y')dx'dy' 
27tr^ 


(3.16) 


in  which  (x,y)  is  the  point  at  which  the  velocity  is  calculated  and  (x',y')  is  an  arbitrary’ 
point  at  which  the  vorticity  is  ^(x',y'). 

The  radial  distance  r  is  given  by: 


r^  =  (X  -  X')-  +  (y  -  y')-  . 


(3.17) 


Normalizing  Eqs.  3.15  -  3.17  through  the  use  of  the  characteristic  dimensions 
given  in  Eq.  3.12 


27ir  * 
m 


(3. IS) 


F  u  =  r 

V  m  J 


27ir  - 

m 


(3.19) 


and 
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(3.20) 


It  is  appropriate  to  take  V^,  the  initial  mutual  induction  velocity  of  the  vortex 
pair,  as  the  characteristic  velocity  and  the  b^,  the  initial  separation  of  the  vorticies.  as 
the  characteristic  length.  Thus,  one  has: 


L 


c 


V 

Fv  =  ^ 


X  b 

C  0 


V 

o 


In  dimensionless  forms,  one  has: 


(3.21) 


N  b 


O  


=  27tF. 


n  =  F  ^  ° 


c  o 


V 


t^  =  tN  =  — ^ 
m  b  F 

o  V 


(3.22) 


in  which  N^b^/V^  is  called  the  stratification  parameter  SP.  The  value  of  SP  normally 
varies  from  zero  to  about  100.  The  larger  the  SP,  the  stronger  the  stratification.  For 
example,  for  a  model  running  at  a  speed  of  5  ft  'sec  at  an  angle  of  attack  (a)  of  10°,  a 
SP‘ smaller  than  1.0  represents  weak  stratification,  a  SP  between  1  and  5  represents 
medium  stratification,  and  a  SP  larger  than  5  corresponds  to  strong  stratification. 

B.  COMPUTER  AIDED  MODELING  OF  VORTEX  MOTION 

The  computer  algorithm  for  the  numerical  solution  of  Eqs.  3.13  and  3.14  is 
provided  in  Apperidix  A.  The  solution  is  done  assuming  the  case  of  inviscid  flow  and 
the  initial  vorticity  distribution  was  assumed  to  be  Gaussian,  i.e.: 


F  '  r  ^  • 
-^exp(-  :^) 


(3.23) 


A  computational  domain  of  5b^  by  lOb^  is  used  in  the  algorithm  in  order  to 
mimick  the  physical  features  of  the  tow  basin.  The’  no-slip  condition  is  satisfied  on  the 
walls  corresponding  to  the  chanel  walls.  The  zero  normal  velocity  condition  is  assumed 
at  the  free  surface  and  along  the  vertical  line  of  flow  symmetiy.  All  computations  are 
done  in  the  fourth  quadrant.  The  fields  of  velocity,  vorticity.  density,  and  the 
fluctuating  components  of  the  density  are  calculated  at  each  time  step  and  plotted  at 
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regular  intervals.  The  computer  algorithm  is  quite  robust  and  can  be  applied  to  any 
stratified  and  unstratified  flow  situation.  As  the  velocity  field  is  computed  from 
analytic  expressions  (i.e.,  the'  Biot-Savart  equations)  vice  an  algorithm  involving  an 
iterative  scheme,  stability  of  the  program  is  enhanced  at  a  slight  expense  in  execution 
time.  An  upwind  finite  differencing  technique  [Ref  20]  was  used  with  an  original  grid 
dimension  of  24  by  44  nodes  with  a  spacing  of0.2ob^.  This  grid  size  was  subsequently 
increased  to  34  by  64  and  44  by  84  in  order  to  examine  the  dependency  of  the  results 
on  the  grid  size. 

The  mesh  coordinates  (x,y)  covering  the  computational  domain  (of  dimension 
m-4  by  n-4)  are: 


X  e  2,3,4,...,m-2 

(3.24) 

Ye  {2,3,4,.. .,n-2}. 

The  necessity  for  the  the  overlap  of  two  nodes  on  all  sides  of  the  computational 
domain  is  to  facilitate  satisfying  the  boundary  condiditons.  Eqs.  3.13  and  3.14  are 
solved  over  the  original  mesh  which  results  in  vorticity  and  density  information. 

The  two  complex  velocity  grids  are  derived  by  shifting  the  original  mesh  h/2  in 
the  X  direction  for  the  u  velocity  and  h/2  in  the  y  direction  for  the  v  velocity.  The 
values  of  these  new  nodes  are  determined  -using  the  vorticity  derived'  for  the  original 
computational  grid  using  the  Biot-Savart  law,  Eqs.  3.15  and  3.16.  The  velocities  at  the 
shifted  grid  points  are  then  averaged  across  the  computational  grid  points  for  further 
use  in  the  finite  differencing  scheme. 
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IV.  DISCUSSION  OF  RESULTS 


The  computer  code  described  in  the  previous  section  has  been  used  to  carry  out 
calculations  for  four  specific  cases.  The  first  of  these  dealt  with  the  case  of  the  rise  of 
vortices  in  a  homogeneous  unstratified  medium.  The  grid  size  as  well  as  the  core 
radius  were  varied  within  limits  to  e.xplore  the  differences  in  the  motion  of  the  vortices. 
Figures  B.l  through  B.13  show  the  velocity  and  vorticity  contours  at  various 
normalized  limes.  Clearly,  the  initial  vonex  spacing  increases  with  time  and  the 
vortices  remain  symmetrical  due  to  the  conditions  imposed  on  the  numerical  analysis. 
In  other  words,  these  vortices  cannot  be  subjected  to  sinusoidal  instability  or  core 
bursting.  Thus,  the  numerical  results  should  agree  more  closely  with  the  experimental 
cases  where  the  vortices  did  not  undergo  such  instabilities.  Figure  B.14  shows  the 
e.xperimental  and  numerical  results.  The  experimental  data  are  separated  as  'HIGH' 
and  'LOW'  because  of  the  fact  that  the  crests  and  troughs  of  the  sinusoidal  instability 
were  separately  evaluated.  The  numerical  results  for  three  normalized  core  sizes  are 
also  shown  in  Figure  B.14.  It  is  clear  that  the  core  size  has  a  pronounced  effect  on  the 
numerical  predictions  and  that  the  results  obtained  with  r^  =  0.08  agree  more  closely 
with  the  experimental  data.  As  noted  earlier,  the  numerical  res-ults  should  be  compared 
with  the  average  instantaneous  height  rather  than  with  the  'HIGH'  and  'LOW'  of  the 
vortices.  There  is  no  conclusive  analysis  or  data  in  the  literature  regarding  the  vortex 
core  radius.  It  has  been  shown  previously  [Ref  16]  that  the  core  radius  depends 
strongly  on  the  shape  of  the  trailing  edges  of  the  lifting  surface.  The  better  rounded 
the  edges  of  the  control  surface  the  tighter  the  winding  of  the  vortex  sheets  , and  thus 
the  smaller  the  vortex  core.  It  appears  that  a  vortex  core  of  0.07  or  0.08  should  be 
used  in  future  calculations. 

Figures  B.15  through  B.31  show  lines  of  constant  density,  velocity,  density 
perturbations,  and  vorticity  contours  for  the  statification  parameter  of  1.0.  F^.  = 
0.02976,  r^  =  0.09,  Y;b„  =  -8.0  and  X.b  =  0.5.  The  srid  used  was  44  bv  84. 
These  figures  show  that  at  small  times  the  vortices  rise  vertically  and  the  circulation  in 
the  flow  field  is  primarily  due  to  the  initial  vortices.  As  time  increases  (see  e.g.,  Figure 
B.21b),  the  regions  of  circulation  with  countersigned  vorticity  develope  in  the  upper 
right  and  lower  left  regions  of  the  vortex.  It  is  easy  to  show  that  this  counter  vorticity 
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not  only  reduces  the  rise  velocity  of  the  original  vortex  pair  but  pushed  the  pair  against 
each  other.  Consequently,  the  countersigned  vorticity  begins  to  dominate  the  flow  (see 
e.g.,  Figure  B.2Sb)  and  the  vortex  pair  migration  stops.  With  further  increases  in  time 
the  vortex  pair  begins  to  migrate  downward. 

The  density  contours  reveal  the  same  phenomenon  in  a  difTerent  context.  .A,s  the 
vortices  rise,  fluid  of  greater  density  is  pushed  upwards  (see  e.g.,  Figure  B.25a)  into 
regions  of  lesser  density.  Since  such  a  migration  cannot  go  on  indefinitely,  the  vortices 
rise  to  a  maximum  height  and  then  begin  to  sink  downwards.  The  calculations  do  not 
take  into  account  the  sinusoidal  instability  and  the  vortex  breakdown.  Consequently, 
the  vortices  in  the  numerical  calculations  continue  to  exist  unimpeded  by  the  instability 
mechanisms.  In  reality,  the  vortices  begin  to  break  up  as  they  near  the  end  of  their 
maximum  migration  and  eventually  disappear. 

The  experimental  and  calculated  values  are  compared  in  Figure  B.32.  The 
correspondence  between  the  measured  and  calculated  values  is  surprisingly  good  up  to 
the  time  of  maximum  rise.  This  is  partly  because  of  the  experimental  fact  that  the  rise 
of  vortices  in  a  highly  stratified  medium  (N  b  /V  =  1.00,  F„  =  0.02976)  is  not 
strongly  dominated  by  sinusoidal  instability  or  vortex  breakdown.  The  migration  of 
vortices  is  inhibited  primarily  due  to  the  reduction  of  vorticity  of  the  initial  vortices  and 
the  creation  of  countersigned  vorticity. 

The  numerical  calculations  were  also  applied  to  a  medium  with  a  relatively 
srnaller  stratification  for  which  the  rise  of  vortices  are  closer  .to  that  of  the 
homogeneous  case  and  thus  the  effect  of  the  core  radius  is  more  pronounced.  The 
results  obtained  with  N'  b^'V  =  0.50,  F„  =  0.02976  and  with  r  =  O.OS  are  shown  in 
Figures  B.33  through  B.43.  The  observations  made  earlier  regarding  the  velocity, 
vorticity,  and  density  contours  may  be  made  with  this  particular  stratification  also. 
Obviously,  the  vortices  rise  to  a  greater  height  and  the  countersigned  vorticity  is 
relatively  weaker  and  developes  at  a  later  stage.  The  calculations  with  this  particular 
stratification  have  been  repeated  by  varying  only  the  initial  core  radius.  Figure  B.44 
shows  the  representative  experimental  data  and  the  numerical  results  obtained  with 
initial  core- radii  of  0.07,  O.OS.  and  0.09.  Evidently  and  as  anticipated  on  the  basis  of 
the  previous  calculations  the  initial  core  radius  of  0.07  yields  results  which  are  in 
excellent  agreement  with  those  obtained  experimentaly.  Figure  B.44  also  shows  that 
the  rise  and  demise  of  vortices  must  be  a  strong  function  of  their  initial  core  radii. 
Experiments  in  both  the  laboratory  and  the  field  have  shown  that  to  be  true  and 
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pointed  out  the  fact  that  the  life  of  a  vortex  is  dictated  partly  by  the  conditions 
prevailing  at  its  creation  and  partly  by  the  conditions  surrounding  the  vortex  at  later 
stages  of  its  motion. 

The  results  presented  so  far  have  dealt  with  vortex  pairs  released  at  relatively 
large  depths  of  submergence  and  have  shown  that  the  computer  code  predicts  results 
which  are  in  good  agreement  with  those  obtained  experimentally  provided  that  a  core 
radius  of  about  r  =  0.07  is  used  in  the  calculations.  It  has  further  been  shown  that 

O 

the  smaller  the  stratification  parameter  the  more  dependent  the  results  become  on  the 
initial  core  readius.  Thus,  the  numerical  analysis  has  not  only  provided  a  power  of 
prediction  for  the  motion  of  the  vortices  but  also  helped  to  delineate  the  proper  value 
of  the  initial  core  radius. 

Having  assertained  the  characteristics  of  the  model  and  the  features  of  the 
velocity,  vorticity,  and  density  variations  it  was  deemed  necessary'  to  explore  problems 
of  greater  practicle  and  operational  concern.  It  is  because  of  this  reason  that  the  depth 
at  which  the  vortices  are  generated  was  chosen  to  be  Y/b^  =  —  2.0  and  the 

subsequent  motion  of  the  vortices  was  explored.  Figures  B.45  through  B.52  show'  the 
densitv,  velocitv,  densitv  fluctuations  and  vorticitv  contours  at  suitable  time  intervals. 
As  the  vortices  approach  the  free  surface  and  begin  to  move  sidew-^ays  under  the 
influence  of  their  images,  the  constant  density  contours  begin  to  intersect  with  the  free 
surface  (see  e.g..  Figure  B.49a)  and  move  sideways  also.  The  intersection  of  these 
density  contours  is  expected  to  give  rise  to  intemal-wave-generated  surface  scars.  The 
prediction  of  the  characteristics  of- these  scars  in  terms  of  the  vortices  is  of  importance 
and  requires  that  the  free  surface  condition  be  changed.  This  change  is  to  allow'  for 
free  surface  deformation  under  the  influence  of  the  tangential  as  well  as  normal 
velocities  induced  at  the  free  surface  by  the  vortices. 

Additional  calculations  are  underway  to  examine  the  effects  of  free  surface 
deformation,  viscous  diffusion  and  ambient  turbulence  for  various  degrees-  of  initial 
stratification.  The  results  presented  herein  and  the  initial  comparisons  are  extremely 
encouraging  and  are  expected  to  lead  to  a  better  understanding  of  this  extremely 
complex  and  challenging  phenomenon. 


V.  CONCLUSIONS 


The  investigation  reported  herein  warranted  the  following  conclusions: 

(1)  The  motion  of  trailing  vortices  generated  by  submerged  bodies  in  a  stratified 
or  unstratified  medium  can  be  effectively  modeled  by  the  developed  computer 
algorithm: 

(2)  The  computer  model  has  aided  in  the  delineation  of  the  proper  value  for  the 
initial  vortex  core  radius; 

(3)  The  initial  vortex  core  radius  and  vorticity  strength  are  the  primary  factors 
governing  the  stability  of  the  numerical  results; 

(4)  The  migration  of  vortices  is  inhibited  primarily  due  to  the  reduction  of  the 
vorticity  of  the  initial  vortices  and  the  creation  of  countersigned  vorticity; 

(5)  The  computer  model  could  be  made  more  realistic  through  the  inclusion  of  a 
deformable  free  surface  in  the  computations. 
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APPENDIX  A 

VORTEX  MOTION  COMPUTER  MODEL 
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PURPOSE : 

THIS  PROGRAM  NUMERICALLY  SIMULATES  THE  RISE  OF  A  TRAILING 
VORTEX  PAIR  SHED  OFF  A  LIFTING  SURFACE  SUBMERGED  TN  A 
DENSITY  STRATIFIED  MEDIUM. 


TO  RUN  THIS  PROGRAM: 

1)  GO  TO  "INPUT  PARAMETERS"  SECTION  OF  CODE  AND  CHANGE  ANY  PROBLEM 
PARAMETERS  (  SP ,  FV,  DT,  ZTMAX,  RO ,  MP ,  NP ,  ISAVE)  AS  DESIRED. 

.2)  SET  MAXIMUM  TIME  IN  DATA  STATEMENT  MAXSTP ; 

EX:  DATA  MAXSTP  /N/ 

(  WHERE  N=  MAXIMUM  NUMBER  OF  TIME  STEPS) 

3)  SET  INTERVALS  FOR  PLOTTING  BY  ENTRIES  IN  DATA  STATEMENT  NPLOT ; 
EX:  DATA  NPLOT/  1,2,8,12,20,40,50,75/ 

(  THIS  PLOTS  OUTPUT  AT  THE  1ST , 2ND , 8TH , 12TH , 40TH , 50TH , AND 
75TH  PROGRAM  TIME  STEP) 

4)  CHOOSE  WHAT  OUTPUT  DEVICE  TO  USE,  EITHER  TEK  618  OR  COMPRES 
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(METAFILE  OUTPUT)  BY  GOING’  TO  GRAPHIC  DEVICE  SECTION  OF  PROGRAI-I 
AND  COMMENTING  OUT  WHICH  OF  TWO  DEVICES  YOU  DONT  WANT  TO  USE. 

5)  COMPILE  PROGRAI^  UNDER  VS  FORTRAN  USING  THIS  COMMAND; 

FORTVS  VORTEX  (OPT (3) 

(  NOTE:  THIS  OPTIMIZES  CODE  AT  LEVEL  3  ) 

6)  GO  TO  APPROPRIATE  TERMINAL .  AI'ID  RUN  UNDER  DISSPLA. 

A.  WHEN  THE  DISSPLA  EXEC  PROMPTS  YOU  ASK  FOR  5  CYLINDERS  OF 
TEMPORARY  DISK  STORAGE  AND  ALSO  YOU  MAY  WISH  TO  ISSUE 

B.  WHEN  THE  EXEC  PROMPTS  YOU  YOU  MAY  WISH  TO  ISSUE: 

FILEDEF  06  DISK  VORIOUT  LISTING  (PERM 
IN  ORDER  TO  REROUTE  THE  PRINTED  OUTPUT  TO  A  DISK  FILE 
NAMED:  VORIOUT  LISTING  A1  . 

C.  PROGRAM  WILL  NOW  RUN  SENDING  OUTPUT  TO  YOUR  DESIGNATED 
DEVICES. 


ADDITIONAL  PROGRAM  NOTES: 

1)  THIS  PROGRAM  COMPUTES  DENSITY  AND  VORTICITY  IN  A  (M-2  X  N-2)  GRID 

BY  SOLVING  THE  BOUYANTLY  SCALED  FORM  OF  THE  NAVIER- STOKES 
EQUATIONS  VIA  AN  "UPWIND  DIFFERENCING"  FINITE  DIFFERENCE  METHOD. 

2)  THE  VELOCITY  FIELD  IS  COMPUTED  USING  A  FORM  OF  THE  BIOT-SAVART 

LAW  AT  EACH  NODAL  POINT  INCLUDING  THE  CONTRIBUTION  OF  THE 
"IMAGE  VORTICIES  "  ADDED  TO  INCLUDE  THE  FREE  SURFACE  EFFECT. 

3)  THIS  VERSION  OF  THE  PROGRAM  INCLUDES  A  DISSPLA  ROUTINE 

TO  PLOT  THE  DENSITY  PERTUBATIONS  , CONSTANT  DENSITY, 

•  VORTICITY,  AND  VELOCITY  FIELD  AT  EACH  SELECTED  TIME  STEP. 

THESE  SELECTIONS  FOR  PLOTTING  ARE  INDICATED  BY  ENTRIES 
IN  THE  ARRAY  NSTE? .  MAXIMUM  NUMBER  OF  TIME  STEPS  IS  INDICATED 
BY  AM  ENTRY  IN  THE  ARRAY  MAXSTEP. 
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THE  COMPUTATIONAL  GRID  (M-2)X(N-2)  IS  USED  FOR  ALL  COMPUTATIONS 
EXCEPT  FOR  VELOCITY  CALLS  TO  THE  3IOT-SAVART  SUBROUTINES. 

TWO  ADDITIONAL  GRIDS,  ONE  FOR  U  AMD  ONE  FOR  V,  ARE  SHIFTED  H/2 
UNITS  IN  THE  X  AND  Y  DIRECTION  RESPECTIVELY.  THE  SUBROUTINES 
COMPUTE  U  AND  V  VELOCITIES  ON  THESE  GRIDS  USING  VORTICITY  FROM 
COMPUTATIONAL  GRID.  THE  VELOCITIES  ARE  THEM  AVERAGED  ACROSS  THE 
COMPUTATION.AL  GRID  POINTS  TO  BE  USED  FOR  THE  FINITE  DIFFERENCING. 

THIS  VERSION  OF  THIS  PROGRAM  WILL  BE  WRITTEN  IN  SINGLE  PRECISION 
AND  IN  AN  OPTIMIZED  FORM  OF  THE  ORIGINAL  CODE. 

THIS  PROGRAM  HAS  THE  OPTION  OF  SAVING  THE  RHO  AND  ZETA  ARRAYS 
AND  OTHER  PROGRAI-I  CONSTANTS  IN  ORDER  TO  RESTART  THIS  SIMULATION 
AT  THE  TIME  A  PREVIOUS  RUN  HAD  FINISHED.  PROGRAM  VORRS  IS 
SPECIFICALLY  WRITTEN  TO  INPUT  THE  DATAFILE  WRITTEN  BY  THIS 
PROGRAM  AND  RESTART  THE  SIMULATION  WITH  THE  SAME  PARAMETERS. 

TO  ACCESS  THIS  SAVE  FEATURE  SET  VARIABLE  ISAVE  =1  .  A  FILEDEF 
OF  FORM  . . .  FILEDEF  08  DISK  <FN>  <FT>  <FM>  . . .  MUST  BE  ISSUED 
AT  THE  APPROPRIATE  TIME  IN  EXECUTION. 

(NOTE:  TO  DISABLE  THIS  SAVE  FEATURE  SET  ISAVE  =0  ) 


TO  CHANGE  THE  MESH  DENSITY  OF  THE  COMPUTATIONAL.  DOMAIN  THE 
USER  SHOULD  DO  THE  FOLLOWING: 

A.  ADJUST  M  AND  M 

B.  ADJUST  INITIAL  VORTEX  POSITION:  MP ,  MP 

C.  ADJUST  THE  DIMENSIONS  OF  THE  W  ARRAY  IN  SUBROUTINES 

CPLOT,  ZPLOT,  AND  DRPLOT  TO  BE  (MW, MW) 

D.  IF  THE  GRID  IS  LARGER  THAN  50  X  100  THEN  ADJUST  THE 

DIMENSION  STATEMENTS  IN  THE  MAIN  PROGRAM  AND  IN 
SUBROUTINES  ZPL,  VPLOT  AMD  V2PLOT  TO  THE  NEW  GRID 
DIMENSION. 
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c 
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Q'k:k'k:k'k7^:k7^:k'k'k:k7^'k:ki^7^:k'k-k:k'k'k7^'^:k7^'k'k'k'k7:7::k’k:k7^:k'k'k'k'kr:'k'k'k7:'k7<::k:k:k'k'k'k'k-^:k7z:ki:'k-k:k 

c 

C  MAJOR  VARIABLES  USED  IN  THIS  PROGRAM; 

C 

C  REAL  VARIABLES; 


C  T  = 

NONDIMENS I ONAL  TIME 

C  DT  = 

TIME  STEP 

C  H  = 

NONDIMENS lONAL  GRID  LENGTH  PARAMETER=  DX  =  DY 

C  RO  = 

NONDIMENS lONAL  VORTEX  CORE  RADIUS 

C  XWIDTH= 

NONDIMENS lONAL  WIDTH  OF  COMPUTATIONAL  AREA 

C  ZTMAX  = 

MAXIMUM  VORTEX  STRENGTH 

C  SP  = 

NONDIMENS lONAL  STRATIFICATION  PARAMETER 

C  FV  = 

NONDIMENS lONAL  FROUDE  NUMBER 

C  BVND  = 

NONDIMENSIONAL  BRUNT-VALSA  FREQUENCY 

C 


C  INTEGER 

VARIABLES; 

o 

11 

NUMBER  OF  NODES  IN  X  DIRECTION 

C  N  = 

NUMBER  OF  NODES  IN  Y  DIRECTION 

o 

11 

X  NODAL  POINT  OF  VORTEX  START  POSITION 

o 

11 

Y  MODAL  POINT  OF  VORTEX  START  POSITION. 

n 

It 

X  ARRAY  DIMENSION 

n 

i 

It 

Y  ARRAY  DIMENSION 

C  NFIG  = 

PLOT  NUMBER 

C  NSTEP  = 

ITERATION  NUMBER 

C 

\ 

C  ARRAY  VARIABLES; 


C  RHO  = 

DENSITY  PERTUBATIONS  AT  EACH  NODE  DUE  TO  DYNAMICS 

C  RHOI  = 

INITAL  VALUES  OF  DENSITY  AT  EACH  NODE  DUE  TO  GRADIENT 

C  RHOT  = 

SUM  OF  BOTH  DENSITY  EFFECTS  AT  EACH  NODE 

C  RHNEW  = 

NEW  VALUE  OF  RHO  AT  NEXT  TIME  STEP 

C  ZETA  = 

VORTICITY  AT  EACH  NODE 

C  ZTNEW  = 

NEW  VALUE  OF  VORTICITY  AT  NEXT  TIME  STEP 
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C  U  =  X  COMPONENT  OF  VELOCITY  AT  EACH  NODAL  POINT 

C  V  =  Y  COMPONENT  OF  VELOCITY  AT  EACH*  NODAL  POINT 

C  U2  =  U  VELOCITYS  COMPUTED  AT  BIOT-SAVART  NODAL  POINTS 

C  V2  =  V  VELOCITYS  COMPUTED  AT  BIOT-SAVART  NODAL  POINTS 

C  X  =  X  COORDINATES  AT  EACH  MODAL  POINT 

C  Y  =  Y  COORDINATES  AT  EACH  MODAL  POINT 

C  XU  =  X  COORDINATE  AT  EACH  U  3.  S.  MODAL  POINT 

C  YU  =  Y  COORDINATE  AT  EACH  U  3.  S.  MODAL  POINT 

C  XV  =  X  COORDINATE  AT  EACH  V  B.  S.  NODAL  POINT 

C  YV  =  y' COORDINATE  AT  EACH  V  B.  S.  NODAL  POINT 

C 

C . . . . . 

c 

C  SUBROUTINES  USED  IN  THIS  PROGRAM: 

C 

C  SUBROUTINE  COUT  -  OUTPUTS  NODE  COORDINATES 

C 

C  SUBROUTINE  OUTPT  -  OUTPUTS  VARIOUS  PROGRAM  ARRAYS  WHEN  CALLED 

C 

C  SUBROUTINE  CPLOT  -  CREATES  ARRAY  OF  DENSITY  AT  EACH  POINT 

C  AND  CALLS  DISSPLA  CONTOURING  ROUTINE  TO  PLOT. 

C 

C  SUBROUTINE  ZPLOT 

C 

c  .  ■  - 

C  SUBROUTINE  DRPLOT 

C 

c 
c 

C  SUBROUTINE  ZPL  -  OUTPUTS  VORTICITY  CONTOURS  AS  A  PRINTER  PLOT. 

C 

C  SUBROUTINE  VPLOT  -  CREATES  COMPLEX  VELOCITY  VECTORS  AT  EACH 

C  NODAL  POINT  AMD  THEN  CALLS  THE  FOLLOWING 

C  SUBROUTINES  TO  PLOT  THEM  ON  DISSPLA: 

C  A)  V2PLOT 


-  CREATES  ARRAY  OF  VORTICITY  AT  EACH  POINT  AND 
CALLS  DISSPLA  COUTOURIMG  ROUTINE  TO  PLOT. 

-  CREATES  ARRAY  OF  DENSITY  PERTUBATIONS  AT  EACH 
POINT  AND  CALLS  DISSPLA  CONTOURING  ROUTINE 

TO  PLOT. 
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C  B)  AKIN 

C  C )  AKAX 

C 

C  SUBROUTINE  COM  -  CALLED  BY  2PL  TO  CONTOUR  VORTICITY  ON  PRINTER. 

C 

c 

C  FUNCTIONS  USED  IN  THIS  PROGRAM: 

C 

C  FUNCTION  VBS  -  COMPUTES  THE  V  COMPONENT  OF  VELOCITY  USING 

C  THE  BIOT-SAVART  LAW. 

C 

C  FUNCTION  UBS  -  COMPUTES  THE  U  COMPONENT  OF  VELOCITY  USING 
C  THE  BIOT-  SAVART  LAW. 

C 

C  FUNCTION  PNEW  -  TAKES  THE  FINITE  DIFFERENCE  TO  STEP  THE 
C  EQUATIONS  IN  TIME. 

C 

DIMENSION  RHOI(50,100) , U( 50 , 100 ) , V(50 , 100 ) , ZETA( 50 , 100 ) 

*  ,2TNEW(50,100) ,RHNEW(50 , 100) ,NPLOT(3) ,X(50) ,Y(100) ,RHO(50 , 100) 

/ 

*  ,XU(50) ,YU(100) ,xy(50) ,YV(100) ,U2(50,100) ,V2(50,100) 

CHARACTER'S  LABEL 

COMMON  WK( 15000) 

COMMON  /PN/  DT,I, J,UF,UFA,UB,UBA,VF,VFA,VB,VBA,VPT,FV 
COMMON  /PARK/  H;MvN',PI  '  •  ;  - 

DATA  MAXSTP/40  / 

DATA  NPLOT/1 , 10,20,30,40,50,60,70/ 

C 


C  . (  CALL  GRAPHIC  OUTPUT  DEVICE  ) 

CALL  COMPRS 
C  CALL  TEK613 

C . 


c 

MM=  50 


MM=  100 


c 


C .  INPUT  PROGRAI^  PARAMETERS  . 

C 

DT=  2.000E0 
T=0.0 

XWIDTH  =5.0 
R0=  0.09E0 
SP=  1.20 
FV=0. 029755 
ZTMIN=  l.OE-8 
M=  44 
N=  84 
ISAVE=  1 

C .  DESIGNATE  NODAL  POINTS  TO  PLACE  INITAL  VORTEX  AT  (MP,NP) 

MP=6 

NP=66 

C - COMPUTE  ADDITIONAL  PARAMETERS - - 

C 

C.....  DENSITY  GRADIENT  IN  Y  DIRECTION  INITALLY  =  DRHO 
DRHOI=-SP*SP*FV*FV 

C .  NONDIMENS I ONAL  B.  V.  FREQUENCY  =  BVND 

BVND=  SP*SP*FV*FV 

C .  MAXIMUM  VORTICITY  AT  ANY  POINT  =  ZTMAX 

ZTMAX=-(FV/(R0^R0)) 

C .  GEOMETRIC  CONSTANT  PI 

PI=  4.0E0*ATAN(1.0E0) 

C .  GRID  LENGTH  PARAMETER  =  H 

H=  XWIDTH/ (M- 4) 

C 

C .  COMPUTE  MISC  M  AND  N  PARAMETERS 

MM1=  M-1 
MP1=  M+1 
MM2=  M-2 
NP1=M+1 
NM1=N-1 
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NM2=  N-2 
MW  =  (2’^M)  -  7 
NW  =  NM2  -  1 
C 

C - IMITALIZE  COORDINATE  SYSTEM  ORIGINS  - 

C 

X(l)=  -H 
Y(l)=  H 
XU(1)=  -l.SEO’^H 
YU(1)=  H 
XV(1)=  -H 
YV(1)=  1.5E0*H 

C 

C  .  COMPUTE  COORDINATES  . . 

C  (  NOTE:  COMPUTATIONAL  GRID  IS  IN  4TH  QUADRANT  +X  -Y  ) 
C 

DO  10  1=  2,MP1 

XU(I)=  XU(I-1)+H 
XV(I)=XV(I-1)+  H 

10  X(I)  =  X(I-l)  +  H 
DO  11  J=  2,NP1 

YU(J)  =  YU(J-1)  -  H 
YV(J)  =  YV(J-l)  -  H 

11  Y(J)  =  Y(J-l)  -  H 
.YHT  =  ABS(Y(N-2)) 

C .  OUTPUT  INITAL  PARAMETERS  OF  THIS  RUN  OF  PROGRAM 

C 

WRITE(6, 1275) 

WRITE (6, 1500) 

WRITE(6,1535)  DT,H 
WRITE(6,1540)  XWIDTH,YHT 
WRITE(6,1550)  SP,FV 
WRITE(6,1560)  M,N 
WRITE (6, 1570)  BVND 
C 
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C .  OUTPUT  COORDINATES  OF  NODES . . . 

WRITE(6,1580)  X(MP),Y(NP) 

WRITE(6,1590) 

CALL  COUT(X,Y,MM,NN) 

C 

C  .  SPECIFY  INITAL  VORTICITY  AND  DENSITY  FIELDS  . 

C 

C  (  FOR  GRID  SYSTEM  X,Y  COMPUTATIONAL  GRID  IN  THE  4TH  QUADRANT  ) 
C  (  TERMS  ARGl  -  ARG4  ARE  IN  QUADRANTS  1  -  '4  ) 


ARGMN=  30.0 
XCENT=  X(MP) 

YCENT=  Y(NP) 

DEM=  2.0E0^R0*R0 
ZETA(MP,NP)=  ZTMAX 
DO  30  J=1,N 

DO  30  1=1, M 
TRM1=0.0E0 
TRM2=0.0E0 
TRM3=0.0E0 
TRM4=0.0E0 

ARG1=  ((X(I)-XCENT)*’^2  +  (Y(  J)+YCENT)**2)/DEN 
IF  (ARGl  .GT.  ARGMN)  GO  TO  12 
TRM1=  EXP (-ARGl) 

12  ARG2=  ( (X(r)+XCENT)’^*2  +  (Y(  J)+YCENT)’^’^2)/DEN 
IF(ARG2  .GT.  ARGMN)  GO  TO  14 
TRI^2=  EXP(-ARG2) 

14  ARG3=  ( (X(I)+XCEMT)**2  +  (Y( J) -YCENT) **2 ) /DEN 
IF(ARG3  .GT.  ARGMN)  GO  TO  16 
TRM3=  EXP(-ARG3) 

16  ARG4=  ( (X(I)-XCENT)’^^2  +  ( Y( J) -YCENT) ^*2 ) /DEN 
IF(ARG4  .GT.  ARGMN)  GO  TO  13 
TRM4=  EXP(-ARG4) 

18  IF((J.EQ.NP)  .AND.  (I  .EQ.MP))  GO  TO  20 

ZETA(I , J)=ZTMAX*(-TRM1  +TRM2  -TRM3  +TRM4) 
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IF(ABS(ZETA(I,J))  .LT.  ZTMIN)  ZETA(I , J)=0 .OEO 
20  RHOI(I,J)=  DRHOI*Y(J) 

RHO(I,J)=  O.OEO 
30  CONTINUE 

WRITE (6, 1250) 

C  CALL  OUTPT(ZETA, 'ZTA(I,  J)  '  ,ra-I,NN,0) 

C  CALL  OUTPT(RHOI, 'ROI(I, J) ' ,MM,NM,0) 

C 

C . '  IMITALIZE  TIME  STEP  COUNTERS . 

NFIG=0 

NSTEP=0 

C 

C*x***x****  ENTRY  POINT  OF  TIME  LOOP 

60  CONTINUE 

C 

C . .  COMPUTE  U  VELOCITY  FOR  ALL  NODES  IN  U  VELOCITY  GRID . 

C 

DO  70  1=  1,MP1 
DO  70  J=  1,N 

70  U2(I,J)=  UBS(XU(I) ,YU(J) ,MM,NN,X,Y,ZETA) 

C 

C . . .  COMPUTE  V  VELOCITY  FOR  ALL  NODES  IN  V  VELOCITY  GRID . 

C 

DO  100  1=  1,M 
DO  100  J=  1,NP1 

100  V2(I,J)=  VBS(XV(I) ,YV(J) ,MM,NN,X,Y,ZETA) 

C 

C .  INTERPOLATE  VELOCITIES  FROM  VELOCITY  GRIDS  TO  NODAL 

C  POINTS  ON  COMPUTATIONAL  GRID. 

VMIN=  l.OE-08 
DO  120  1=  1,M 
DO  120  J=  1,N 

U(I,J)=  (U2(I,J)  +  U2(I+1, J))/2.0E0 
V(I,J)=  (V2(I,J)+  V2(I, J+1))/2.0E0 

IF  (ABS(U(I,J))  .LT.  VMIN)  U(I,J)=  O.OEO 
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IF  (ABS(V(I,J))  .LT.  VMIN)  V(I,J)=  O.OEO 
120  CONTINUE 

C  IF  (NSTEP  .EQ.  0)  CALL  OUTPT (V, ' V( I , J)  ',MM,NN,1) 

C  IF  (NSTEP  .EO.  0)  CALL  OUTPT(U, ' U(I , J)  ',MM,NN,1) 

C 

C 

C  . . 

C  INCREASE  TINE  COUNTER  3Y  ONE. 

C  CHECK  ITERATION  COUNTER  AND  TERMINATE  IF  NSTEP  IS  MORE  THAN  MAX 
C  COMPUTE  VORTICITY  AND  DENSITY  AT  ALL  INTERIOR  POINTS 

C  (  STORED  AS  ELEMENTS  OF  ZTNEW  AND  RHNEW  ) 

C . 

c 

NSTEP=NSTEP+1 

IF  (NSTEP  .GT.  MAXSTP)  GO  TO  1700 
T=  T+DT 

DO  600  1=2, MMl 

DO  600  J=2,NM1 
VPT=  V(I,J) 

UF=  (U(I+1,J)+U(I, J))/2.0E0 
UFA=ABS(UF) 

UB=  (U(I-1, J)+U(I,J))/2.0E0 
UBA=ABS(UB) 

VF=  (V(I, J-1)+V(I, J))/2.0E0 
VFA=ABS(VF) 

VB=  (V(I,J+1)+V(I,J))/2.0E0 
VBA=ABS(VB) 

ZTNEW ( I , J ) =  PNEW ( ZETA , RHO , MM , NN , - 1 . OEO , 0 . OEO , 0 . OEO ) 


600  RHNEW(I,J)=  PNEW(RHO,RHO,MM,MN,O.OEO,O.OEO,BVMD) 

C 

C 

C - APPLY  BOUNDRY  CONDITIONS  TO_  RHO  AND  ZETA . . 

C 

C .  ASSIGN  NEW  VALUES  OF  RHO  AMD  ZETA  TO  RHO  AND  ZETA  ARRAYS  . . 

C 
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c 


DO  700  I  =  2,MM1 
DO  700  J  =  2,NH1 

2ETA(I,J)=  ZTMEW(I,J) 

700  RHO(I, J)=RHMEW(I , J) 

C .  CALCULATE  RHO  AND  ZETA  VALUES  AT  I=M-2  AMD  J=N-2  BOUNDRIES  ... 

C 

DO  300  1=  2,HM2 

DO  300  J=  2, MHZ 

ZETA(MH2,J)=  (ZETA(MM1,J)+  ZETA(M-3, J))/2.0E0 
ZETA(I,NH2)=  (ZETA(I,NM1)+  ZETA( I ,N-3 ) ) /2 . OEO 
RHO(MM2,J)=  (RH0(MH1,J)+  RHO (M-3 , J) ) /2 . OEO 
RHO(I,NH2)=  (RH0(I,NH1)+  RHO ( I ,N-3 ) )/2 . OEO 

IF(ABS(ZETA(I,J)) .LT.  ZTMIN)  ZETA(I , J)=0 .OEO 
IF(ABS(RHO(I, J) ) .LT.  ZTMIN)  RHO(I,J)=  O.OEO 
BOO  CONTINUE 
C 

C . 

C  PLOT  FLOW  PATTERN  WHEN  NSTEP  =  VALUE  SPECIFIED  IN  NPLOT 
C 

C  THEN  RESTART  ALGORITHM  AT  LINE  60  FOR  A  NEW  TIME  STEP 

C.  . . . . . . . . . 

DO  1100  K=  1,3 

IF  (NSTEP-NPLOT(K) )  1100,  1200,  1100 

1100  CONTINUE 
GO  TO  60 
1200  NFIG=NFIG+1 
TM=  T’^FV 

IF  (MSTE?  .EQ.  HAXSTP  )  CALL  OUTPT(ZETA, ' ZTA(I ,J) ', MM, NN, NSTEP) 
IF  (NSTEP  .EQ.  MAXSTP)  CALL  OUTPT  (RHO,  • '  RHO  ( I ,  J)  ,  MM  ,NN  ,NSTEP ) 
WRITE (6, 1300 )NSTEP,T,TM 

CALL  CPLOT ( RHO , RHOI , MM , NM , NSTEP , TM , MW , MW ) 

CALL  VPLOT ( U , V , X , Y , MM , NN , TM ) 

CALL  DRPLOT ( RHO , MM , NM , NSTEP , TM , MW , NW ) 

CALL  ZPL(ZETA,MM,NN,XWIDTH,YHT, ITER) 
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CALL  2PL0T(ZETA,MI-I,NN,TM,MW,NW) 

1250  FORMAT(//  10X,'INITAL  ZETA  AND  RHO  VALUES  :  ') 

1275  FORMAT( ' 1 ' ) 

1300  FORI-IAT(/10X,  '  ITERATION  MUM  : ',16 , 4X ,  '  TIME= '  ,  F9 . 3 , 4X ,  '  T’^  =’,F10.6) 
1400  WRITE (6, 1500) 

1500  FORMAT (  ///) 

1530  FORMATC  ',5X,'INITAL  PROGRAM  PARAMETERS:'  ) 

1535  FORMAT( 'O' ,5X, 'TIME  STEP=  ' , ?3 . 4 , 5X , ' GRID  SPACING  H  =  ',”8.4) 
1540  FORMAT( 'O' ,5X, 'NON  DIMENSIONAL  CELL  DIMENSIONS  :  XDIRECTION  =', 

*  F8.4,5X, 'YDIRECTION  =',F8.4) 

1550  FORMAT('0' ,5X, 'STRATIFICATION  PARAMETER  =',F8.4,5X, 

*  'FROUDE  NUMBER  =',F8.4) 

1560  FORMAT( 'O' ,5X, 'NUMBER  OF  NODES  IN  X  DIRECTION  =',I4,5X, 

*  ' IN  Y  DIRECTION  =' ,14) 

1570  FORMAT( 'O' ,5X, 'NONDIMENSIONAL  3.  V.  FREQUENCY  =',F15.9) 

1580  FORMAT( 'O' ,5X, 'X  COORDINATE  OF  VORTEX  START  POSITION  =',F8.4,5X, 
'Y  COORDINATE  =',F8.4) 

1590  FORMAT( 'O' ,5X, 'OUTPUT  COORDINATES  OF  ALL  OTHER  NODES:') 

C 

C  LOOP  BACK  TO  STARTING  POINT  OF  ITERATION 

GO  TO  60 

1700  CALL  DONEPL 
C1700  CONTINUE 
C 

C .  SAVE  ARRAYS  AND  TIME  DATA  TO  RESTART  PROGRAM  LATER  . 

C 

IF  (ISAVE  .EQ.  0  )  GO  TO  1900 
WRITE(3,*)  SP,FV,DT,T 
WRITE ( 8 , ^ )  M , N , MP , NP 
WRITE (8,^)  BVND,H,ZTMAX 
DO  1800  1=  1,M 

DO  1800  J=  1,N 

1800  WRITE(8,*)  ZETA(I,J),  RH0(I,J) 

1900  STOP 
END 
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c 


REAL  FUNCTION  PNEW (P , Q, MM , NN , A , B , C) 

C 

C  THIS  FUNCTION  INTEGRATES  THIS  EQUATION  WITH  RESPECT  TO  TIME 
C  $$$  DP/DT=  (-D(UP)/DX  -D(VP)/DY  +  A*DQ/DX  +  B’^DEL(P)  +C*V  )  $$$ 

C  UTILIZING  AN  UPWIND  DIFFERENCING  SCHEME  ' 

C 

DIMENSION  P(HM,ra),Q(MM,NN) 

COMl^ON  /PN/  DT,I, J,UF,UFA,UB,UBA,VF,VFA,VB,VBA,VPT,FV 
COra-ION  /PARM/  H,M,N,PI 

Pl=  (UF-UFA)*P(I+1,J)  +(UF+UFA-UB+UBA)’^P(I,  J)  -  (UB+UBA)’^P(I-1 ,  J) 
P2=  (VF-VFA)*P(I, J-1)  +(VF+VFA-VB+VBA)*P(I, J)  - (VB+VBA)*P ( I , J+1 ) 
P3=  Q(I+1,J)  -Q(I-1,J)  . 

P4f  P(I+1,J)  +P(I-1,J)+P(I,J-1)+P(I,J+1)-4.0E0’^P(I, J) 

P5=  2.0E0*H*VPT 

TEMP1=(DT/(2.0E0*H))*((-P1-P2)  +A*P3+ (2 . OEO/H) *B*P4+C*P5 ) 
PNEW=P(I, J)+TEMP1 

•  RETURN  ...  .  •  ■ 

END 

C 

C==================================================================== 

c 

REAL  FUNCTION ' UBS ( XU , YU , MM , NN , X , Y , ZETA ) 

C 

C  THIS  SUBROUTINE  FINDS  FOUNDRY  VALUES  OF  VELOCITY  U 
C  UTILIZING  THE  BIOT-SAVART  LAW 

C 

C  (NOTE:  TERMS  IN  EQUATION  ARE  EVALUATED  CONSECUTIVELY  IN  EACH 


41 


c 


QUADRANT  FROM  ONE  TO  FOUR  ) 


C . 

DIMENSION  X(MM) , Y(NN) , ZETA(MM,NN) 

COMMON  /FARM/  K,M,N,PI 
C 

UHIN=  l.OE-08 
Ul=  O.OEO 
U2=  O.OEO 
U3=  O.OEO 
U4=  O.OEO 
C 

SIGNP=1 .OEO 
SIGNM=-1 .OEO 
DA=  H*H 
C 

C  . -  COMPUTE  U  COMPONENT  OF  VELOCITY  USING  THE  B-S  LAW - 

c  -  ■ 

DO  20  1=  1,M 

DO  20  J=  1,N 
C  --  QUADRANT  ONE  -- 

RSQ1=((XU-X(I))**2  +  (YU+Y(J))**2) 

Ul=  U1+  ((-Y(J)-YU)/(2.0E0*PI*RSQ1))*SIGNM*2ETA(I,J)*DA 
C  --  QUADRANT  TWO  -- 

RSQ2=  ((XU+X(I))**2  +(YU+Y(J))*’^2) 

U2=  U2+  ((-Y(J)-YU)/(2.0E0^PI*RSQ2)  )’^SIGNP*2ETA(I,  J)*DA 
C  --  QUADRANT  THREE  -- 

RSQ3=  ((XU+X(I)  )*’^2  +(YU-Y(J)  )**2) 

U3=  U3+((Y(J)-YU)/(2,0E0’^PI’^RSQ3)  )’^SIGNM^ZETA(I,  J)*DA 
C  --  QUADRANT  FOUR  -- 

RSQ4=  ((XU-X(I))’^^2  +  (YU-Y(J)  )’^*2) 

20  U4  =  U4  +  ((Y(J)-YU)/(2.0E0*PI*RSQ4))*SIGMP*ZETA(I,J)’^DA 
UTEMP=  UH-U2  +U3+  U4 
IF  (ABS(UTEHP) .LT.UHIN)  UTEMP=0.0E0 
UBS=  UTEMP 
RETURN 
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END 


C 

C 

c================================================================== 

c 

REAL  FUNCTION  VBS (XV , YV ,MM ,NN , X, Y , ZETA) 

C 

C.  THIS  SUBROUTINE  FINDS  BOUMDRY  VALUES  OF  VELOCITY  V 
C  UTILIZING  THE  BIOT-SAVART  LAW 

C 

C  (NOTE  :  TERNS  IN  THE  EQUATION  ARE  EVALUATED  CONSECUTIVELY  IN  EACH 
C  QUADRANT  FROM  ONE  TO  FOUR) 

C . 

DIMENSION  X(MM) ,Y(NN) ,ZETA(MN,NN) 

COMMON/PARM/  H,M,N,PI 
C 

VMIN=:  l.OE-08 
Vl=  O.OEO 
V2=  O.OEO 
V3=  O.OEO 
.  V4=  O.-OEO 

C  '  .  ■ 

SIGNP=1.0E0 
SIGNN=-1 .OEO 

.DA=  H*H  •  • 

C 

c  -  COMPUTE  V  COMPONENT  OF  VELOCITY  USING  THE  B-S  LAW - 

C 

DO  20  1=  1,M 

DO  20  J=  1,M 
C  --  QUADRANT  ONE  -- 

RSQ1=((XV-X(I))**2  +  (YV+Y(J))**2) 

Vl=  V1+  ((XV-X(I))/(2.E0*PI^RSQ1))*SIGNT'I*ZETA(I,J)*DA 
C  --  QUADRANT  TWO  -- 

RSQ2=  ( (XV+X(I)  )**2  +(YV+Y(  J)  )*’^2) 


43 


V2=  V2+  ((XV+X(I))/(2.0E0’^PI*RSQ2))*SIGNP’^ZETA(I,J)’^DA 
C.  --  QUADRANT  THREE  -- 

RSQ3=  ((XV+X(I))’^’^2  +(YV-Y(J))’^*2) 

V3=  V3+  { (XV+X ( I ) )  /  ( 2 . 0E0*PI’^RSO3 ) )  *SIGMM’^ZETA(  I ,  J )  ’^DA 
C  --  QUADRANT  FOUR  -- 

RSQ4=  ( (XV-X(I)  )’^’^2  +  (YV-Y(J))^’^2) 

20  V4  =  V4  +  ((XV-X(I))/(2.0E0*PI*RSQ4))*SIGNP’^ZETA(I,J)*DA 
VTEMP  =  V1+V2  +V3+  V4 
IF(ABS(VTSMP) .LT.VNIN)  VTENP=0.0E0 
VBS=  VTENP 
RETURN 
END 
C 

c 

C - - - - - - - - - 

c 

SUBROUTINE  COUT(X,Y,MN,NM) 

C 

c 

- - 

REAL  X(MM) ,Y(NN) ,H,PI 
COMMON/PARM/H , M , N , PI 
WRITE(6,50) 

WRITE (6,*)  '  ARRAY  X' 

■  ••  WRITE(6,*) 

WRITE(6,*)  (X(I),  1=1, M) 

WRITE(6,*) 

WRITE(6,*) 

WRITE (6,*)  '  ARRAY  Y 

WRITE  (6,*) 

WRITE(6,*)  (Y(J),  J=1,N) 

WRITE(6,^) 

WRITE (6,*) 

50  FORMAT (///) 

RETURN 
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END 


C===================================================================: 

SUBROUTINE  OUTPT ( A , LBL , MM , NN , KAL ) 

C 

C  THIS  SUBROUTINE  PRINTS  OUT  AN  (NKM)  ARRAY  WITH  A  LABEL 
C  IN  A  RECTANGULAR  (  2 , M- 2 )X( 2 ,N-2 )  GRID 

C 

DIMENSION  A(MM,MN) 

COMMON/PARM/  K , M / N , P I 
CHARACTER’^8  LBL  ' 

NM2=  N-2 
MM2=  M-2 
WRITE (6, 50) 

WRITE (6,*) 

WRITE (6,*)  '  OUTPUT  OF  ARRAY:  ',LBL,'  ITER  NUM KAL 

WRITE (6,*) 

DO  10  J=  2,NM2 
WRITE  (6,=^) 

WRITE(6,*)  (A(I,J),  1=2, MM2) 

10  CONTINUE 
50  FORMAT( '1' ) 

RETURN 

END 

C 

C- . . . . . - . . . . . . - . . . 

c 

SUBROUTINE  CPLOT ( RHO , RHOI , MM , NN , NSTEP , TM , MW , NW) 

COMMON  WK( 15000) 

COMMON  /FARM/  H,M,N,PI 
REAL  W(81,81) 

REAL  RHO(MM,NN) ,RHOI(MM,NM) , H , PI , RHOT( 50 , 100) 

DATA  ISCALE/4HSCAL/ 

C 

C... CREATE  ARRAY  OF  DENSITY  CHANGES  AT  EACH  NODE  (DRHO)  . 

C 
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MM2=  M-2 
MM2=  N-2 
DO  10  1=  2,HM2 

DO  10  J=  2,NH2 
JJ=  (NM2+1)-J 
11=  I+(M-5) 

RHOT ( I , J ) =RHO ( I , J ) +RHO I ( I , J ) 

w(ii,jj)=  r:-:ot(i,j) 

10  CONTINUE 

DO  20  1=  3, MM2 
11=  (M-l)-I 

DO  20  J=  2,NM2 
JJ=  (NM2+1)  -J 


20  W(II,JJ)=  RHOT(I,J) 

C 

C  CALL  OUTPT ( RHOT , ' RHT ( I , J ) J , MM , NN , NSTEP ) 

C .  CALL  SUBROUTINE  CONTR  TO  PLOT  THESE  PERTU3ATI0NS 

C 


CALL  RESET  (3HALL) 

CALL  PAGE(8.7,11.2) 

CALL  PHYSOR(2.0,3.00) 

CALL  AREA2D  (5. 0,6.0) 

CALL  HEIGHT (.2  ) 

CALL  INTAXS 

‘CALL'XNAME  ('X/BO',4) 

CALL  YNAME  ( 'Y/BO' ,4) 

CALL  XTICKS(2) 

CALL  YTICKS(2) 

CALL  GRAF  ( - 5 . 0 , 1 . , 5 . , 10 . , 1 . , 0 . ) 

CALL  MIXALF( ' INSTRU' ) 

CALL  MESSAG('T(EH.5)^(EXHX)  =$',100,1.5,6.5) 
CALL  REALNO(TM, 106 , 'ABUT' , 'ABUT' ) 

CALL  FRAME 

CALL  BCOMON  (15000) 

CALL  CONMIN(2.0) 
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CALL  CONANG  (20.) 

CALL  CONMAK  ( W , MW , MW , ISCALE ) 

CALL  HEIGHT (.08) 

CALL  CONLIN  ( 0 , 5HSOLID  ,  '  LABELS  ' 2 , 5 ) 

CALL  CONLIN  ( 1 , 4HDASH , 3HNOLABELS , 1 , 3 ) 

CALL  RASPLN(0.25) 

CALL  COMTUR  ( 2 , 6HLABELS , 4HDRAW) 

CALL  HEIGHT (.2) 

CALL  COMPLX 

CALL  RESET ( 'HEIGHT' ) 

CALL  RESET ('COMPLX') 

CALL  ENDPL  (0) 

RETURN 
END  ' 

C 

C . - . . . . 

c 

SUBROUTINE  ZPLOT ( ZETA , MM , NN , TM , MW , NW ) 

REAL  ZETA(MM,NN) ,H,PI 
COMMON  WK( 15000) 

COMMON  /PARM/  H,M,N,PI  .  , 

DIMENSION  W(81,81) 

DATA  ISCALE/4HSCAL/ 

C 

C... CREATE  ARRAY  FOR  PLOTTING  OF  ZETA  IN  ACTUAL  COMPUTATIONAL  AREA.. 
C 

m2=  M-2 
NM2=N-2 

DO  10  1=  2, MM2 
11=  I+(M-5) 

DO  10  J=  2,NM2 
JJ=  (NM2+1)-J 
10  W(II,JJ)=ZETA(I,J) 

DO  20  1=  3, MM2 
11=  (M-l)-I 
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c 

c 

c 

c 


DO  20  J=  2,NM'2 
JJ=  (m2+l)  -  J 
W(II,JJ)  =  -ZETA(I,J) 

CALL  SUBROUTINE  COMTR  TO  PLOT  VORTICITY 


CALL  RESET  (SHALL) 

CALL  PAGE (8. 7, 11. 2) 

CALL  PKYSOR(2.0,3.00) 

CALL  AREA2D  (5. 0,6.0) 

CALL  HEIGHT (.2  ) 

CALL  INTAXS 
CALL  XMAME  ( ' X/BO ' ,4) 

CALL  YNAllE  (  ‘Y/BO  '  ,4) 

CALL  XTICKS(2) 

CALL  YTICKS(2) 

CALL  GRAF  ( -5 . 0 , 1 . , 5 . , 10 . , 1 . , 0 . ) 

CALL  MIXALF( ‘INSTRU' ) 

CALL  MESSAG( 'T(EH.5)’^(EXHX)  =$',100,1.5,6,5) 
CALL  REALNO(TM,106, 'ABUT' , 'ABUT' ) 

CALL  FRAME 

CALL  BCOMOM  (15000) 

CALL  CONAMG  (60.) 

CALL  CONMIN(2.0) 

CALL  CONMAK  (W,MW,NW, ISCALE) 

CALL  HEIGHT (.03) 

CALL  CONLIN  (0 , 5HSOLID , ' LABELS ' , 2 , 5 ) 

CALL  CONLIN  ( 1 , 4HD ASH ,' LABELS ', 1 , 3 )  ' 

CALL  RASPLN(0.25) 

CALL  CONTUR  ( 2 , 6HLABELS , 4HDRAW) 

CALL  HEIGHT ( .2) 

CALL  COMPLX 

CALL  RESET ( 'HEIGHT' ) 

CALL  RESET ( ' COMPLX' ) 
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CALL  ENDPL  (0) 

RETURN 

END 

C 

C 

C-- . - . . 

c 

SUBROUTINE  DRPLOT ( RKO , MM , MM , MSTE? , TM , MW , MW ) 

COMMON  WK( 15000) 

COmON  /FARM/  H,M,N,PI 
REAL  W(81,81) 

REAL  RHO(MM,NN) ,H,PI 
DATA  ISCALE/4HSCAL/ 

C 

C... CREATE  ARRAY  OF  DENSITY  CHANGES  AT  EACH  NODE  (DRHO) 
C 

MM2=  M-2 
NM2=  N-2 
DO  10  1=  2, MM2 

DO  10  J=  2,NM2 
JJ=  (NM2+1)-J 
11=  I+(M-5) 

W(II,JJ)=  RHO(I,J) 

10  '  CONTINUE 

DO  20  I=, 3, MM2 
■  11=  (M-l)-I 

DO  20  J=  2, MM2 
JJ=  (NM2+1)  -J 


20  W(II,JJ)=  RHO(I,J) 

C 

C .  CALL  SUBROUTINE  CONTR  TO  PLOT  THESE  PERTUBATIONS 

C 


CALL  RESET  (SHALL) 
CALL  PAGE(8.7,11,2) 
CALL  PHYSOR(2.0,3.00) 
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CALL  AREA2D  (5. 0,6.0) 

■  CALL  HEIGHT( .2  ) 

CALL  IMTAXS 

CALL  XNAME  ('X/BO' ,4) 

CALL  YMAME  ( ' Y/BO ' ,4) 

CALL  XTICKS(2) 

CALL  YTICKS(2) 

CALL  .GRAF  (-5.0,1.,5.,10.,1.,0.) 

CALL  MIXALF( ' INSTRU' ) 

CALL  MESSAG('T(EH.5)*(EXHX)  =$',100,1.5,6.5) 

CALL  REALNO (TM , 106 , ' ABUT ' , ' ABUT ' ) 

CALL  FRAME 

CALL  BCOMOM  (15000)  • 

CALL  CONMIN(2.0) 

CALL  CONAMG  (20.) 

CALL  CONMAK  (W,MW,^^■/,  ISCALE) 

CALL  HEIGHT (.08) 

CALL  CONLIN  (0 , 5HSOLID , ' LABELS ' , 2 , 5 ) 

CALL  CONLIN  ( 1 , 4HDASH , 8HNOLABELS , 1 , 3 ) 

CALL  RASPLN(0.25) 

CALL  CONTUR  ( 2 , 6HLABELS , 4HDRAW) 

CALL  HEIGHT (.2) 

CALL  COMP LX 

CALL  RESET ('HEIGHT') 

CALL  RESET ( 'COMPLX' ) 

CALL  ENDPL  (0) 

RETURN 

END 

C 

0====================================================================== 

c 

SUBROUTINE  CON ( 22 , XRANGE , YRAMGE , NX , H , MM , NN , ZMAX , ZMIN ) 

C 

C  THIS  SUBROUTINE  GENERATES  A  CONTOUR  PLOT  OF  ZETA  IN  A  RECTANGULAR 

C  DOMAIN  OF  SIZE:  YRANGE  *  XRANGE  ON  THE  PRINTER. 
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c 


CHARACTER’^  1  SYMBOL (8)  , GRAPH (i 00) 

REAL  LX,LY 

DIMENSION  Z2(MM,NN) , VALUE (7) 

DATA  SYl-IBOL  /  '  X '  ,  '  A '  ,  '  -  '  ,  '  '  ,  '  '  ,  '  +  7 

VALUE (1)=  ZMIN 

VALUE(2)=  ZMIN*. 7 5 

VALUE (3 )=  ZMIN*. 5 

VALUE(4)=  0.0 

VALUE (6 )=  ZMAX*.75 

VALUE (7 )=  ZMAX 

VALUE(5)=ZMAX*.45 

NY=  0 . 8*NX* YRANGE / XRANGE 

DELX=  XRANGE/NX 

DELY=  YRANGE/NY 

WRITE (6, 50) 

WRITE ( 6 , * )  ' NX= ' , NX ,  ' NY= ' , NY , ' XRANGE= ' , XRANGE , 

*  '  YRANGE=  '  ,  YRAI'IGE 

WRITE (6,*) 

DO  11  1=  1,NX 
11  GRAPH(I)=  'T' 

WRITE(6,6)  (GRAPH(I),  1=  1,NX) 

DO  5  JSYMBL  =  1,NY 
Y=  (JSYMBL-0.5)*DELY 
J=  1+Y/H 
LY=  Y-(J-1)*H 
HMLY=  H-LY 

DO  4  ISYl^BL  =  1,NX 

X=  (ISYM3L-0.5)*DELX 

I=l+X/H 

LX=  X-(I-1)*H 

HMLX=  H-  LX 

Al=  HMLX*HMLY 

A2=  HMLX*LY 

A3=  LX*LY 


51 


A4=  LX’^HMLY 

Z2C=  (A1*Z2(I,J)  +A2’^Z2(I,  J+1)  +A3*Z2(I+1,  J+1) 

*  +A4=*'Z2(I+1,  J))/H*’^2 
C 

C...  DETERMINE  THE  VALUE  OF  MKAMGE  BASED  OM  Z2C . 

C 

C -  PRINT  OUT  CHARACTERS  IN  ARRAY  GRAPH  BY  ROW  . 

C 

DO  2  K=  1,7 

IF  (Z2C-VALUE(K))  1,  1,  2 

1  NRANGE  =  K 

GO  TO  3 

2  CONTINUE 
NRANGE =  8 

3  GRAPH (  ISYMBL)=  SYMBOL (NRANGE ) 

4  CONTINUE 

5  WRITE(6,6)  (GRAPH(I),  1=1, NX) 

6  FORMAT ('  ',10X,100A1  ) 

50  FORMAT ( ' 1 ' ) 

DO  60  1=  1, NX 
60  GRAPH(I)=  'B' 

.  WRITE(6,6)  (GRAPH(I)  ,1  =  1, NX)  . 

WRITE(6,*) 

WRITE (6,*)  'MAXIMUM  VALUE  OF  NEGATINE  VORTICITY  =  SYMBOL  (X)  = 

*  ,ZMIN  .  .  - 

WRITERS,*)  '  75%  OF  MAXIMUN  NEGATIVE  VORTICITY  =  SYMBOL  (A)  ' 
WRITE (6,^)  '  50%  OF  MAXIMUM  NEGATIVE  VORTICITY  =  SYMBOL  (-)  ' 
WRITE (6,*) 

WRITE (6,*)  '  MAXIMUM  VALUE  OF  POSITIVE  VORTICITY  =  SYMBOL  (1)  = 
,  ZMAX 

WRITE(6,^)  '  75%  OF  MAXIMUM  POSITIVE  VORTICITY  =  SYMBOL  (+) ' 

RETURN 

END 

C . . . . 

c 
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SUBROUTINE  ZPL (ZETA , MM , NN , XWIDTH , YHT , ITER ) 

COMMON  /FARM/  H,M,N,PI 
DIMENSION  Z2(50,100) ,ZETA(MM,MN) 

C 

C... CREATE  ARRAY  FOR  PLOTTING  OF  ZETA  IN  ACTUAL  COMPUTATIONAL  AREA.. 
C 

MN3=  M-3  . 

MM3=N-3 
ZMIN=  0.0 
ZMAX=  0.0 
DO  10  1=  1,MM3 

DO  10  J=  -1  ,NM3 
Z2(I,J)=  ZETA(I+1,J+1) 

C..  FIND  THE  MAXIMUM  POSITIVE  VORTICITY  IN  ARRAY  Z2(MM,NN) . 

C 

IF(Z2(I,J)  .GT.  ZMAX)  ZMAX=Z2(I,J) 

C..  FIND  THE  MAXIMUM  NEGATIVE  VORTICITY  IN  ARRAY  Z2(MM,NN) . 

C 

IF(  Z2(I,J)  .LT.  ZMIN)  ZMIN=  Z2(I,J) 

10  ■  CONTINUE 

C 

C  ‘  • 

IF  (ZMAX  .LT.  l.OE-6  )  ZMAX=1.0 

C .  CALL  SUBROUTINE  CON  TO  PLOT  VORTICITY . 

c  ■ 

CALL  CON ( Z2 , XWIDTH , YHT , 1 0 0 , H , MM , NN , ZMAX , ZMIN ) 

RETURN 

END 

C- . 

SUBROUTINE  VPLOT ( U , V , X , Y , MM , NN , TM ) 

COMMON  /FARM/  H,M,N,PI 
COMPLEX  ZV( 50, 100) ,  CV(50,100) 

DIMENSION  U(MM,NM),  V(MM,NN),  X(MM) ,  Y(NN) 

MM2=  M-2 
NM2=  N-2 
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MM3  =  M  -  3 
NM3  =  N  -  3 
C 

C .  FORM  COMPLEX  U,  V,  AND  Z  ARRAYS  FOR  QUADRANTS  3  AND  4 

C 

DO  10  1=  2,MH2 

DO  10  J=  2,NM2 
11=  I-l 

JJ=  (NH2+1)  -J 
ZV(II,JJ)=  CMPLX(X(I) ,Y(J) ) 

10  CV(II,JJ)=  CMPLX(U(I,J),V(I,J)) 

C 

C 

C .  PLOT  VECTORS  USING  SUBROUTINE  V2PLOT  . 

CALL  V2PLOT(ZV,CV,MM3,NM3,TM) 

C 


RETURN 

END 


C 

C . 

C  • 

SUBROUTINE  V2PLOT (Z , C , M ,N ,TM) 


C 

C 

C 

C 

C 

C 


A  X  ^  A  X  X  X  X  A  A  X  x  A  7C  X  A  X  A  A  X  X  X  x  X  X  X  A  A  A  X  7^  x  X  7^  x  7^  X  X  x  X  t:  X  X 

*  SUBROUTINE  V2PLOT  TO  PLOT  MXN  COMPLX  VELOCIT  FIELD  C(M,N)  ^ 

*  OVER  A  COMPLX  DOMAIN  Z(M,N)  * 

t:  a 


7VA7*:A:^A7*CAAA7^AA7*C7*C7*CA7^7*r7:AAAAAAA7V:*C7^AAAAA:*CAA7^A:^A7^AAXAA7^A7*:7^A7^7^:*CA7*CAAA7^ 


COMPLEX  Z(50,100) ,C(50,100) ,CMAX 
REAL  X(IOO) ,Y(100) ,YX(100) ,YN( 100) ,XX(100) , 
*  XN(IOO) ,X2(100) ,Y2(100) 

INTEGER  PMAX,PMIN 
EXTERNAL  AMAX 
EXTERNAL  AMIN 


C 
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DO  20  J  =1,N 


DO  10  I  =1,M 

X(I)  =ABS(REAL(C(I, J))) 
Y(I)  =ABS(AIMAG(C(I,J))) 
10  CONTINUE 

CALL  AMAX(Y,M,YMAX,PMAX) 
CALL  AMAX(X,M,XI'IAX,?MAX) 
YX(J)  =  YMAX 
XX(J)  =  XMAX 
CALL  AMIN(Y,N,YNIN,PMIN) 
CALL  AMIN(X,M,XMIN,PMIN) 
YN(J)  =YNIN 
XN(J)  =XMIN 
20 _  CONTINUE 

•  CALL  AMAX(YX,N,YMAX,PMAX) 

C  CALL  AMIN(YN,N,YNIN,PMIN) 

CALL  AMAX(XX,N,XMAX,PMAX> 

C  CALL  AMIN(XN,N,XHIN,PHIN) 

C 

CMAX  =  CMPLX( XMAX, YMAX) 

UMAX  =SQRT(  (XMAX)  ^^^^2+ (YMAX)  **2) 
DO  24  I  =1,M 
DO  22  J=1,M 

C(I,J)  =  C(I,J)/UMAX 
22  CONTINUE 

24  CONTINUE 

C 

DO  40  J  =1,N 

DO  30  I  =1,M 

X(I)  =REAL(2(I,J)) 

Y(I)  =AIMAG(2(I,J)) 

30  CONTINUE 

CALL  AMAX ( Y , M , YMAX , PMAX ) 
CALL  AMAX (X,M, XMAX, PMAX) 
YX(J)  =  YMAX 
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XX(J)  =  XMAX 

CALL  AMIN(Y,M,YMIN,PMIN) 

CALL  AMIN(X,M,XMIN,PMIN) 

YW(J)  =Y!-IIM 
XN(J)  =XHIN 
40  CONTINUE 

CALL  AMAX ( YX , N , YMAX , PMAX ) 

CALL  AH I N ( YN , N , YH I N , PHI N ) 

CALL  AHAX(XX,N,XHAX,PHAX) 

CALL  AMIN(XN,N,XHIN,PMIN) 

C 

CALL  RESET (SHALL) 

CALL  PAGE (8. 7, 11. 2) 

CALL  PHYSOR(2.0,3.0) 

CALL  AREA2D(5. ,6.) 

CALL  HEIGHT (.2) 

CALL  INTAXS 

CALL  XNAHE( 'X/BO ' ,4) 

CALL  YNAHE( 'Y/BO' ,4) 

CALL  XTICKS(2) 

CALL  YTICKS(2) 

CALL  GRAF(0 .0,0. 5, 5. 0,-10. 0,1. 0,0.0) 

CALL  HIXALF( ' INSTRU' ) 

CALL  HESSAG( 'T(EH.5)*(E:^.X)  =$’,100,1.5,6.5) 
CALL  REALNO(TM, 106, 'ABUT' , 'ABUT' ) 

C 

DO  60  J=1,N 

DO  50  1=1, M 

X(I)  =  REAL(Z(I,J)) 

X2(I)  =  X(I)+REAL(C(I, J)) 

Y(I)=  AIHAG(Z(I,J)) 

Y2(I)  =  Y(I)+AIHAG(C(I,J)) 

50  CONTINUE 

CALL  HARKER  (3) 

CALL  SCLPIC( .001) 
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54 

60 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


c 


5000 

5100 


CALL  CURVE (X,Y,M,-1) 

DO  54  I  =1,M 
.  X11=X(I)  • 

X12=X2(I) 

Y11=Y(I) 

Y12=Y2(I) 

IF(X12.EO.O. .AND.Y12.EQ.0)  GO  TO  54 
CALL  RLVEC  (Xll ,Y11 ,X12 ,Y12 , 1121) 
CONTINUE 
CONTINUE 

.  CALL  RESET ( 'HEIGHT' ) 

CALL  ENDPL(O) 

RETURN 

END 


SUBROUTINE  AMAX(Y,N, YMAX,PHAX) 

7*:  x  X  7^  7^  X  7C  7^  X  7^  X  7^  7^  7^  7^  X  X  7k  X 

*  SUBROUTINE  TO  COMPUTE  THE  LARGE  ELELEMENT 

*  IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  ' A ' . 

7k  7k  7k  7k  7k  7k  7k  7k  7k  X  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  7k  X  7k  7k  7k 

***  VARIABLE  DECLARATION 
REAL  AMAX  ‘  . 

INTEGER  PMAX 
DIMENSION  Y(N) 

YMAX=Y(1) 

DO  5100  1=2, M 

IF{Y(I) .LE.YMAX)  GO  TO  5000 
YMAX=Y(I) 

PMAX=I 
CONTINUE  • 

CONTINUE 
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c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


6000 

6100 


MAXROW=POS 

RETURN 

END 


SUBROUTINE  ANIN ( Y , N , YMIN , PMIN ) 


X  X  X  7C  X  X  7C  7?  X  X  7?  x  X  X  7^  X  X  X  X  X  X  t:  X  TT  7^:  :r  TT  X  X  x  X  X 


*  SUBROUTINE  TO  COMPUTE  THE  SMALL  ELEMENT 

*  IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  ' A ' . 

7^  7^  7^  7^  7^  7^  7^  ^  7^  ^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  7^  T:  7^  7^  7^  X  X  7^  7^  7^ 


***  VARIABLE  DECLARATION 
REAL  AMIN 

INTEGER  PMIN 
DIMENSION  Y(N) 

YMIN=Y(1) 

DO  6100  1=2, N 

IF(Y(I) .GE.YMIN)  GO  TO  6000 
YMIN=Y(I) 

PMIN=I 

CONTINUE 

CONTINUE 

RETURN 

END 
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Y/BO 


APPENDIX  B 
FIGURES 


Figure  B.la  Velociiy  Field:  SP  =  0.00,  T*  =  0.05952. 
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Y/BO 


0.05952 


I  ^  i  ^  i  ^  ^  ^  ^  ^  ^  '  i  ^  i  '  i  ’  ^  ' 

-5-4-3  -2  -1  0  1  2  3  4  5 


X/BO 

Figure  B.lb  Vorticity  Field;  SP  =  0.00. 
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Y/BO 


T'  =  OY  59520 


Figure  B.2a  Velocity  Field:  SP  =  0.00. 
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Y/BO 


T'*  =  0  =  59520 


-5-4-3  -2  -1  0  1  2  3  4  5 


,X/BO 


Figure  B.2b  Vorticity  Field;  SP  =  0.00. 
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T'  =  1  =  19040 


Figure  B.3a  Velocity  Field:  SP  =  0.00. 
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Y/BO 


1 .19040 


O 


Figure  B.3b  Vorticity  Field:  SP  =  0.00. 
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1  .78560 


Figure  B.4a  Velocity  Field:  SP  =  0.00. 
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Y/BO 


1.78560 


T 


K 
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CD  - 


CO  - 


CD  -i 
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I  ^  ^  ^  ^  ^  i  ^  I 

-3  -2-10  1 

X/BO 


Figure  B.4b  Vorticity  Field;  SP 


1  ‘  I  '  I  ' 

2  3  4  5 

0.00. 
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=  2.38080 


Figure  B.5a  Velocity  Field;  SP  =  0.00. 
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Y/BO 


T’'  =  2.38080 


O 


CM  - 


X/BO 


Figure  B.5b  Vorticity  Field:  SP  =  0.00. 
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Y/BO 


2=97600 


Figure  B.6a  Velocity  Field:  SP  =  0.00. 
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Y/BO 


2.97600 


X/BO 


Figure  B.6b  Vorticity  Field;  SP  =  0.00. 
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T'  =  3.57120 
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Figure  B.7a  Velocity  Field:  SP 
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Figure  B.7b  Vorticity  Field:  SP  =  0:00. 
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Figure  B.8a  Velocity  Field;  SP  =  0.00. 
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Figure  B.8b  Vorticity  Field;  SP  =  0.00. 
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Figure  B.9a  Velocity  Field:  SP  =  0.00. 
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Figure  B.9b  Vorticity  Field;  SP  =  0.00. 
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Figure  B.lOa  Velocity  Field:  SP  =  0.00. 
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Figure  B.lOb  Vorticity  Field:  SP  =  0.00. 
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Figure  B.lla  Velocity  Field:  SP  =  0.00. 


79 


Y/BO 


T’'  =  5.95200 


X/BO 


Figure  B.llb  Vorticity  Field:  SP  =  0.00. 
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Figure  B.12a  Velocity  Field:  SP  =  0.00. 


81 


Y/BO 


1“'  =  6.54720 


X/BO 


Figure  B.l2b  Vorticity  Field;  SP  =  0.00. 
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Figure  B.13a  Velocity  Field:  SP  =  0.00. 
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Figure  B.13b  Vorticity  Field:  SP 
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Figure  B.14  Experimental  and  Numerical  Results  for  Various  Initial  Core  Radii. 
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Fieure  B.15a  Constant  Densitv  Contours:  SP  =  1.00. 
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Figure  B.15b  Velocity  Field:  SP 
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Figure  B.15c  Density  Perturbation  Contours:  SP  =  1. 00. 
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Figure  B.l5d  Vorticity  Field;  SP  =  1.00. 
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Figure  B.16a  Constant  Density  Contours;  SP 


1.00. 


90 


=  0.34380 


X/BO 


Figure  B.16b  Velocity  Field:  SP  =  1.00. 
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Figure  B.16c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.16d  Vorticity  Field:  SP  =  l.OO. 
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Figure  B.17a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.17b  Velocity  Field:  SP  =  1.00. 
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Figure  B.17c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.17d  Vorticity  Field:  SP  =  1.00. 
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Figure  B.lSa  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.lSb  Velocity  Field:  SP  =  1. 00. 
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Fiaure  B.18d  Vorticitv  Field:  SP  =  1.00. 
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Figure  B.l9a  Constant  Density  Contours;  SP  =  1.00. 
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Figure  B.I9b  Velocity  Field:  SP  =  1.00. 
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Figure  B.19c  Density  Perturbation  Contours:  SP  =  1.00. 


104 


Y/BO 


1 .37520 


K 


■5  -4  -3  -2-10  1 

X/BO 


2  3 


5 


Figure  B.19d  Vorticity  Fieid:  SP  =  1.00. 


105 


Y/BO 


O 


T'  =  1.71900 


Figure  B.20a 


Constant  Density  Contours:  SP 


1.00. 


106 


=  1  .71900 


0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5 

X/BO 

Figure  B.20b  Velocity  Field;  SP  =  1.00. 


107 


Y/BO 


1 .71900 


T 


-5-4-3-2-10  1  2  3  4 


5 


X/BO 


Figure  B.20c  Density  Perturbation  Contours;  SP  =  1.00. 
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Figure  B.20d  Vorticity  Field;  SP  =  1.00. 
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Figure  B.21a  Constant  Density  Contours;  SP  =  1.00. 
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Figure  B.21c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.21d  Vorticity  Field:  SP  ==  1.00. 
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Figure  B.22a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.22b  Velocity  Field:  SP  =  1.00. 
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Figure  B.22c  Density  Perturbation  Contours:  SP  =  l.OO. 
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Figure  B.23a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.23b  Velocity  Field:  SP  =  1.00. 
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Figure  B.23c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.23d  Vorticity  Field:  SP  =  1.00. 


121 


Y/BO 


T**  =  3-.09420 


-5  -4  -3  -2  -1  0 


2  3  4  5 


X/BO 


Figure  B.24a  Constant  Density  Contours:  SP  =  1.00. 


122 


r  = -3=09420 


^  ^  ^  ^  ^ 

^  J,  -  /.  /.  '- 

^  .  7  . 

^  ^  ^  ^  ^  ^ 

^  ^  ^  ^ 

^  ^  ^  ^ 
^  ^  ^  ^ 


-^.-C^<<TT•V^^ 

-4  r  r  '^  ■ 

-6.  <  ^  x 


K  A  A  A  A  /I  . 

A  A  A  A  /  /  > 

A  A  A  A  ^  ^ 

A  n  'I  ' 
'I'l  ^  ‘ 
'i''  '}  'I  - 1  A  A  A  ‘ 

A  AAAA'AAAA' 

A  A.  'i  'i  'i  ^  rf  ^  ^  ^  y 

:v  ‘i  ^  5l  1  i  ^  ^  r  r  r  r 

^a!  ^Af  ^  ^  ^  ^  ^  ’T  ^  “T  “T  ^ 


^^^rTT'-r-TT-r 
T^^r-r-rrT'-r 

^  ^  ^  i\  i  A 


\ 

\ 

“ 

\ 

f  1 

CD  J 

1  ! 

\ 

ll 

i  1 

o  i 

'  ' 

! 

r - 1  > 

f  'i 

LJ 

\>  \>  \x  \x  ^ 

\>i^-S^-^>4’7'i  A  K  K  Kv^V^v^^^^VYt-r’TT'r’r'r' 

-1  A  A  K  K  ^^^'f'^^J'^^^■^'^r'r'r■ 
A  A  A  K  K  S 

^  ^  ^  b  r  r  r  ^^^^^T't- 

A  A  A  A  h  K  KK 

^3:5V'S>.^^^^>-^5''^'?’‘?’'^-^7'f4-l4  A  A  h  h  K  ^^^^^ 


X/BO 


Figure  B.24b  Velocity  Field:  SP  =  1.00. 
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Figure  B.24c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.24d  Vorticity  Field:  SP  =  1.00. 
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Figure  B.25a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.25b  .Velocity  Field:  SP  =  1.00. 
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Figure  B.25c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.25d  Voriicity  Field:  SP  =  1.00. 
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Figure  B.26a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.26b  Velocity  Field;  SP  =  l.OO. 
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Figure  B.26c  Density  Perturbation  Contours:  SP  =  1.00. 


132 


Y/BO 


3.78180 


X/BO 


Figure  B.26d  Vorticity  Field:  SP  =  1.00. 
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Figure  B.27a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.27b  Velocity  Field:  SP  =  1.00. 
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Figure  B.27c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.27(l  Vorticity  Field:  SP  =  1.00. 
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Figure  B.2Sa  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.28b  Velocity  Field;  SP  =  1.00. 
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Figure  B.28c  Density  Perturbation  Contours:  SP  =  1.00. 


140 


Y/BO 


4=46940 


-5  -4  -3  -2  -1  0 


2  3  4  5 


X/BO 


Figure  B.28d  Vorticity  Field:  SP  =  1.00. 
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Figure  B.29a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.29b  Velocity  Field:  SP  =  l.OO. 
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Figure  B.29c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.29d  Vorlicity  Field:  SP  =  l.OO. 
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Figure  B.30a  Constant  Density  Contours:  SP  =  1.00. 
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Figure  B.30b  Velocity  Field;  SP  =  1.00. 
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Figure  B.30c  Density  Perturbation  Contours:  SP  =  1.00. 
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Figure  B.31a  Constant  Density  Contours;  SP  =  1.00. 
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Figure  B.SIb  Velocity  Field:  SP  =  1.00. 
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Figure  B. 31c  Density  Perturbation  Contours:  SP  =  1. 00. 
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Figure  B.3ld  Vorticiiy  Field:  SP  =  l.OO. 
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Figure  B.32  .  Experimental  and  Numerical  Results  for  SP  =  1.00. 
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Figure  B.33a  Constant  Densitv  Contours:  SP  =  0.50. 
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Figure  B.33b  Velocity  Field;  SP  =  0.50. 
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Figure  B.33c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.33d  Vorticicy  Field;  SP  =  0.50. 
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Figure  B.34a  Constant  Density  Contours;  SP  =  0.50. 


159' 


Y/BO 


T^=  0.59520 


Figure  B.34b  Velocity  Field;  SP  =  0.50. 
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Figure  B.34c  Density  Perturbation  Contours;  SP  =  0.50. 
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Figure  B.34d  Vorticity  Field;  SP  =  0.50. 
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Figure  B.35a  Constant  Density  Contours:  SP  =  0.50.. 
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Tigure  B.35b  Velocity  Field;  SP  =  0.50. 
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Figure  B.35c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.35d  Vorticity  Field;  SP  =  0.50. 
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Figure  B.36a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.36b  Velocity  Field:  SP  =  0.50. 
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Figure  B.36c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.36d 


Vorticity.  Field;  SP  =  0.50. 
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Figure  B.37a  Constant  Density  Contours:  SP  =  0.50. 


171 


T’'  =  2.38080 


Figure  B.37b  Velocity  Field:  SP  =  0.50. 
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Figure  B.37c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.37d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.3Sa  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.38b  Velocity  Field;  SP  =  0.50. 
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Figure  B.3Sc  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.38d  Vorticily  Field:  SP  —  0.50. 
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Figure  B.39a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.39b  Velocity  Field:  SP  =  0.50. 
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Figure  B.39c  Density  Perturbation  Contours;  SP  =  0.50. 
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Figure  B.39d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.40a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.40b  Velocity  Field:  SP  =  0.50. 
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Figure  B.40c  Density  Perturbation  Contours:  SP 


0.50. 
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Figure  B.40d  Vorricicy  Field:  SP  =  0.50. 
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Figure  B.4la  Constant  Density  Contours:  SP  =  0.50. 


187 


T”  =  4.76160 


Figure  B.41b  Velocity  Field;  SP  =  0.50. 
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Figure  B.41c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.41d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.42a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.42b  Velocity  Field:  SP  =  0.50. 
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Figure  B.42c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.42d  Vorticity  Field;  SP  =  0.50. 
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Figure  B.43a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.43b  Velocity  Fieid;  SP  =  0.50. 
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Figure  B.43c  Density  Perturbation  Contours;  SP  =  0.50. 
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Figure  B.43d  Vorlicity  Field;  SP  =  0.50. 
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Figure  B.44  Experimental  and  Numerical  Results  for  Various  Initial  Core  Radii. 
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I'igure  B.45a  Coiisiant  Density  Contours:  SB  =  0.50. 
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Figure  B.45b  Velocity  Field;  SP  =  0.50. 
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Figure  B.45c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.45d  Voriicity  Field:  SP  =  0.50. 
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Figure 'B. 46a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.46b  Velocity  Field:  SP  =  0.50. 
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Figure  B.46c  Density  Perturbation  Contours;  SP  =  0.50. 
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Figure  B.46d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.47a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.47b  Velocity  F'ield;  SP  =  O.fO. 
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Figure  B.47c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.47d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.4Sa  Constant  Density  Contours;  SP  =  0.50. 


212 


Y/BO 


1.78560 


Figure  B.4Sb  Velocity  Field:  SP  =  0.50. 
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Figure  B.48c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.4Sd  Voriicity  Field:  SP  =  0.50. 
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Figure  B.49a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B,49b  Velocity  Field;  SP  =  0.50. 
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Figure  B.49c  Density  Perturbation  Contours;  SP  =  0'.50. 
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Figure  B.49d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.50a  Constant  Density  Contours;  SP  =  0.50. 
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Figure  B.50b  Velocity  Field:  SP  =  0.50. 
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Figure  BJOc  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.50d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.51a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.51b  Velocity  Field:  SP  =  0.50. 
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Figure  B.51c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.51d  Vorticity  Field:  SP  =  0.50. 
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Figure  B.52a  Constant  Density  Contours:  SP  =  0.50. 
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Figure  B.52b  Velocity  Field:  SP  =  0.50. 
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Figure  B.52c  Density  Perturbation  Contours:  SP  =  0.50. 
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Figure  B.52d  Vorticity  Field;  SP  =  0.50. 
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